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Abstract 

Mechanisms underlying the emergence of orientation selectivity in the primary visual cortex are highly debated. Here 
we study the contribution of inhibition-dominated random recurrent networks to orientation selectivity, and more 
generally to sensory processing. By simulating and analyzing large-scale networks of spiking neurons, we investigate 
tuning amplification and contrast invariance of orientation selectivity in these networks. In particular, we show how 
selective attenuation of the common mode and amplification of the modulation component take place in these 
networks. Selective attenuation of the baseline, which is governed by the exceptional eigenvalue of the connectivity 
matrix, removes the unspecific, redundant signal component and ensures the invariance of selectivity across different 
contrasts. Selective amplification of modulation, which is governed by the operating regime of the network and 
depends on the strength of coupling, amplifies the informative signal component and thus increases the 
signal-to-noise ratio. Here, we perform a mean-field analysis which accounts for this process. 

Keywords: Orientation selectivity; Contrast invariance; Inhibition-dominated; Mean-field analysis; Common-mode 
attenuation; Tuning amplification; Operating regime 



Introduction 

Neurons in sensory cortices of mammals often respond 
selectively to certain features of a stimulus. A well-known 
example in the visual system is the orientation of an elon- 
gated bar (Hubel and Wiesel 1962, 1968). This specificity 
of neuronal responses is believed to be a fundamental 
building block of stimulus processing and perception in 
the mammalian brain. Although well studied for more 
than half a century now, it is not fully clear which neuronal 
mechanisms generate this selectivity. In particular, it is not 
clear which kind of network structure is necessary for its 
establishment, and whether different system architectures 
are being employed by different species. 

In the center of the current debate is the role of feedfor- 
ward vs. recurrent connectivity in the initial establishment 
of selectivity, as well as its further properties like contrast 
invariance (Alitto and Usrey 2004; Niell and Stryker 2008; 
Sclar and Freeman 1982) and tuning sharpening (Ferster 
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and Miller 2000; Sompolinsky and Shapley 1997). Mod- 
els relying on a purely feedforward structure, as it was 
originally suggested by Hubel and Wiesel (1962), cannot 
explain why the orientation tuning of both neuronal spik- 
ing and membrane potentials is invariant with respect to 
stimulus contrast (Anderson et al. 2000) (see Finn et al. 
(2007), however, for a more elaborate feedforward model 
to account for contrast invariance). On the other hand, the 
prevailing recurrent network models for the intra-cortical 
origin of contrast-invariant orientation selectivity (Ben- 
Yishai et al. 1995; Somers et al. 1995) cannot explain how 
highly selective neuronal responses emerge in mice (Niell 
and Stryker 2008), as they rely on feature specific connec- 
tivity which rodents seem to lack, at least at the onset of 
eye opening (Ko et al. 2013) a . Also, it is not clear whether 
the orderly arrangement of preferred features on the cor- 
tical surface, which has been described in most primates 
and carnivores (Blasdel and Salama 1986; Bonhoeffer and 
Grinvald 1991; Ohki et al. 2006; Ts'o et al. 1990), is nec- 
essary for the emergence of feature selectivity, or serves 
any function at all (Horton and Adams 2005). Different 
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answers to these questions would have radically differ- 
ent implications for understanding higher brain functions, 
and how to study them. 

Here we investigate the emergence of contrast invariant 
orientation selectivity in large-scale networks of spiking 
neurons with dominant inhibition. The biological motiva- 
tion for studying such networks come from experimental 
studies which show functional dominance of inhibition 
(Haider et al. 2013; Rudolph et al. 2007). The dominance 
of inhibition is probably a consequence of the dense, 
local pattern of inhibitory connectivity, which has been 
reported for different cortices (Fino and Yuste 2011; Hofer 
et al. 2011; Packer and Yuste 2011). From a theoreti- 
cal point of view, such networks are well-studied (Brunei 
2000; van Vreeswijk and Sompolinsky 1996) and they have 
been shown to exhibit asynchronous -irregular (Al) activ- 
ity states that are in many respects resembling the spiking 
activity recorded in the mammalian neocortex. 

We first show that highly selective tuning curves, which 
are contrast invariant, can be obtained in these networks, 
even in absence of any feature-specific connectivity and 
any spatial network structure. We then analyze these net- 
works by proposing a simplified mean-field description, 
which predicts the main properties of output orienta- 
tion selectivity in the networks. The analysis identifies 
the mechanisms responsible for tuning amplification and 
contrast invariance. We show that the results hold for a 
wide range of parameters, and for networks operating in 
different recurrent regimes. 

Results 

Contrast invariant orientation selectivity in random 
networks 

We consider a large-scale network model of spiking neu- 
rons, representing a small volume of cerebral cortex. Our 
model consists of a recurrent network of 12 500 leaky 
integrate-and-fire (LIF) neurons, 80% excitatory and 20% 
inhibitory (Braitenberg and Schiiz 1998). Each neuron 
receives input from 10% of the excitatory population and 
10% of the inhibitory population, sampled randomly. This 
is the same network configuration as the one considered 
by Brunei (2000). The network is strongly inhibition dom- 
inated, as individual inhibitory synapses are arranged to 
be 8 times more effective than the excitatory ones (see 
Methods for details). We fix the random connectivity, and 
only change the strength of recurrent coupling, measured 
as the amplitude of the postsynaptic potential (PSP) at 
each synapse. The post-synaptic currents are modeled as 
alpha-functions, with time constant r syn (see Methods). 
We refer to the peak value of PSP, J a , as EPSP, and to the 
total PSP as / = T syn e/ a . The default set of parameters is 
given in Table 1. 

The response of the network with EPSP = 0.1 mV to 
the stimulus of certain orientation, 6, is shown in Figure 1. 



Table 1 Table of notations and parameters 



Neuron model 


Membrane time constant 


r m 


20 ms 


Resting potential 


^rest 


OmV 


Threshold voltage 


VtU 


20 mV 


Reset voltage 


Preset 


20 mV 


Refractory period 


Tref 


2 ms 


Synaptic model 


Synaptic time constant 


T syn 


0.5 ms 


Peak EPSP 


EPSP 


0.01 1 mV 


Synaptic efficacy 


J 


r syn e EPSP 


Inhibition dominance ratio 


g 


8 


Feedforward strength 


EPSPfiw 


0.1 mV 


Synaptic delay 


D 


1.5 ms 


Network connectivity 


Number of neurons 


N 


12 500 


Excitatory fraction 


f 


0.8 


Connection probability 


€ 


10% 


Weight matrix 


W 


Wlj 


Network operator 


A 


a - w)-i 


Simulation 


Stimulus orientation 


e 


0°, 15° 165° 


Preferred orientation (PO) 


6* 


[0°,180°) 


Contrast 


c 


9,39,99% 


Baseline firing rate 




12,16, 20 kHz 


Modulation ratio 


m 


m = 1 0% 


Simulation time 


fsim 


6000 ms 


Analysis 


Orientation selectivity index 


OSI 




Orientation selectivity vector 


osv 




Scatter degree index 


SDI 




Baseline (modulation) gain 


Kb (Ym) 




Linearized neuronal gain 


? 




Tuning curve 


W) 





The external input, s, that a neuron receives is modeled 
as a homogeneous Poisson process, which is slightly (m = 
10%) modulated by the orientation of the stimulus accord- 
ing to a cosine function: s(9) = sg [1 + mcos(2(6* — 9*))]. 
A random preferred orientation (PO), 6*, of the external 
input is assigned to each neuron. The feedforward efficacy 
is fixed in all simulations to EPSPff w = 0.1 mV. The raster 
plot of the activity for one second is shown in Figure 1A, 
which indicates that with these parameters the network is 
indeed operating in the asynchronous irregular (Al) state. 
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Figure 1 Raster plot of network activity. (A) Typical raster plot, shown is 1 second of the activity of a balanced network with recurrent synaptic 
couplings of amplitude EPSP = 0.1 mV, in response to a stimulus of orientation 8 = 90° at a medium contrast (sg = 16 000 spikes/s). Spikes of 
excitatory and inhibitory neurons are plotted in red and blue, respectively. The plot on the right shows the average firing rate of neurons, computed 
from the spike count during 6 s of simulation. The histogram on the bottom depticts the time resolved firing rates of the excitatory (red) and the 
inhibitory (blue) populations, respectively, using a time window of 10 ms.The histogram on the bottom right shows the probability density function 
of time averaged firing rates of individual excitatory (red) and inhibitory (blue) neurons in the network, respectively. (B) Same spike trains as above, 
with neurons sorted according to their input preferred orientations (indicated on the y-axis). As above, the plot on the right indicates the firing rate 
of each neuron computed from the spikes emitted during the simulation. 



If neurons are sorted according to their preferred ori- 
entation, the differences in the firing rates become visible 
(Figure IB). Neurons with a preferred orientation closer 
to the orientation of the stimulus on average respond with 
higher rates, while the neurons closer to the orthogonal 
orientation are mostly silent. The cosine tuning of the 
input is reflected by the cosine tuning of firing rates across 
the population (Figure IB, right). 

If we repeat the stimulation of the network with dif- 
ferent orientations, the individual tuning curves for each 
neuron are obtained. For a sample neuron from the same 
network this is shown in Figure 2A-B. Although the neu- 
ron receives input that is only weakly tuned, the network 
is capable of amplifying the selectivity, and the output tun- 
ing is much more pronounced. This emerging selectivity 
is independent of stimulus contrast, C, reflected by the 
stimulus-specific change in the mean firing rate of the 
external input, sb- Figure 2A shows, for a sample neuron, 
that the shape of tuning curves remains unchanged for dif- 
ferent contrasts. Normalizing the output tuning curves by 
the peak value (Anderson et al. 2000) yields exactly the 
same curve for all input intensities (Figure 2A, inset). The 
other neurons in the network show the same behavior, as 
the polar plots in Figure 2C demonstrate for a larger sam- 
ple. Note that excitatory and inhibitory neurons are both 
highly selective. 

To quantify orientation selectivity across the popula- 
tion, we compute two orientation selectivity indexes, OSI 



(Niell and Stryker 2008; Ringach et al. 2002), for all tuning 
curves (Figure 3). Both measures (Figure 3A and B) show 
that the mean OSI is increased by the network, to a level 
that is comparable with the results reported in animal 
experiments (Chapman and Stryker 1993; Dragoi et al. 
2001; Niell and Stryker 2008; Ringach et al. 2002). More- 
over, the selectivity is maintained upon increasing the 
input intensity, as both the mean value and the shape of 
OSI distributions are very similar for different contrasts. 

To directly verify this invariance, we compare the OSI of 
all neurons at different contrasts. Plotting the OSIs at the 
medium contrast (MC) vs. the lowest contrast (LC), and 
at the highest contrast (HC) vs. the medium, indeed reveal 
that the majority of neurons show a remarkable robust- 
ness of their tuning curves upon a change in contrast 
(Figure 3C). 

The high selectivity in the network emerges despite the 
fact that each neuron receives input from a large pool 
of neurons with heterogeneous selectivity and different 
preferred orientations (Figure 3D). In fact, the PO dis- 
tribution of presynaptic neurons is essentially uniform 
(Figure 3D, top histogram), and the presynaptic OSI dis- 
tribution (Figure 3D, right histogram) is very similar to 
the OSI distribution of the whole population (Figure 3A). 
Therefore, the output response is highly selective, despite 
the fact that the input is quite heterogeneous, as reported 
in experiments (Chen et al. 2011; Jia et al. 2010; Varga 
et al. 2011). 
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Figure 2 Tuning amplification and contrast invariance of neuronal responses in the network. (A) The random network (center), composed of 
excitatory (red, 80%) and inhibitory (blue, 20%) neurons, transforms weakly tuned input (green) to highly selective responses (red). The sample 
tuning curves are shown for unit # 4715, which is an excitatory neuron. Neurons have very similar tuning curves for different contrasts, and they are 
indistinguishable after normalization by their respective peak value (see inset). EPSP = 0.1 mV. (B) Spikes elicited by the same excitatory neuron as 
in (A) in response to different orientations of a stimulus at a medium contrast. (C) More samples of output tuning curves from the network, in polar 
representation. Radial axis encodes the firing rate of the respective neuron at each orientation, indicated by the angle, with the maximum firing rate 
indicated as a number next to each plot. Both excitatory (red) and inhibitory (blue) neurons are highly selective in their responses, despite their 
weakly selective inputs (green) and recurrent inputs of mixed selectivities. Lighter colors correspond to tuning curves for lower contrasts, respectively. 



This result is similar to a recent study of orientation 
selectivity in rodents (Hansel and van Vreeswijk 2012), 
in that both show random networks are capable of gen- 
erating selective output responses. In the following, we 
provide a detailed mathematical analysis of the mecha- 
nisms involved in this process. Our mean-field analysis 
indeed enables us to compute the mean output responses 
of networks quite precisely. 



single-neuron gain. If v = 0, the input-output relation is 
given by 



-f+Wr + J s s, 



(2) 



and if the matrix 1 



r=L(t 



W)~ l ~s = As. 



W is invertible, this readily implies 

(3) 



A reduced linear rate model of the network 

To illustrate the main network processing, we first recruit 
a linear rate model of the network. We start by a reduced 
diagrammatic description of the network (Douglas et al. 
1995). This is equivalent to the description of the network 
in its stationary state, in terms of neuronal firing rates r 
that arise as a result of a stimulus 5 driving a recurrent 
network W (see Methods, Section Mean firing rates in 
recurrent networks): 

v=-T+W? + J s s. (1) 

The matrix W encodes the recurrent synaptic connec- 
tions in a network comprising N neurons, s and r are the 
AT-dimensional vectors of firing rates of input/stimulus 
and output/response, respectively, v is the AT-dimensional 
vector of time-averaged membrane potentials, and J s is the 



The single-neuron gain J s in Eq. (2) subsumes the feedfor- 
ward processing of input to all neurons, before recruiting 
any lateral interactions (a in Figure 4A). The operation 
of the recurrent network on a rate vector r is given by 
the feedback Wr appearing in the same equation. If the 
feedback gain /J was the same for all activity configura- 
tions, y — — P) would be the associated closed- 
loop gain of the system (Figure 4A). In this case, any 
stimulus would be amplified or attenuated by the net- 
work with a uniform factor y. As a consequence, if (p 
is a stimulus feature systematically varied in an exper- 
iment, the output tuning curves r(4>) would have the 
same shape as the input tuning curves s((p), as they 
would merely be rescaled by y as a whole. However, the 
amplification of orientation selectivity observed in our 
simulated networks suggests that the amplification fac- 
tor of the untuned part (baseline) is considerably smaller 
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Figure 3 Orientation selectivity across the population. (A) Distribution of a global measure of orientation selectivity, OSI = 1 — Circular Variance 
(Ringach et al. 2002), in the network. Lighter colors show the distributions for lower contrasts, respectively. All inputs have an OSI of 0.05. 
(B) Distribution of an alternative measure of orientation selectivity often used by experimentalists (Niell and Stryker 2008). OSI* is the difference of 
activity at preferred and orthogonal orientations, normalized by their sum, (r pte f — r 0 rth)/(fpref + forth)- (/pref ar| d r orth) are obtained from the best fit 
of a cosine function to output tuning curves, evaluated at Output PO and Output PO + 90°, respectively. Alternatively, OSI* can be computed from a 
linear interpolation of data points (inset). Lighter colors show the distributions for lower contrasts, respectively. All inputs have an OSI* of 0.1. (C) The 
OSI of all neurons for medium contrast (MC) vs. low contrast (LC), and for high contrast (HC) vs. medium contrast are plotted in the left and right 
panels, respectively. The diagonal line indicates a perfect contrast invariance of OSI. (D) Output PO vs. Output OSI for all neurons in the presynaptic 
pool of the neuron shown in Figure 2A. A stimulus of medium contrast has been applied. The neuron receives input from presynaptic neurons that 
are themselves highly selective on average (OSI distribution on the right), and which uniformly cover the whole range of possible Output POs 
(distribution on top). The Output OSI and Output PO of the target neuron are 0.65 and 105°, respectively. Other neurons receive similarly 
heterogeneous inputs (not shown). 



than that of the tuned part (modulation). As a conse- 
quence, the selective amplification/attenuation of the net- 
works considered here is reflected by an extended diagram 
(Figure 4B) with two separate channels, and different over- 
all gains for baseline and modulation, respectively. Even 
if the feedforward gain is identical for baseline and mod- 
ulation, the feedback gain is different for each stimulus 
component, leading to different closed-loop gains yb and 
Ym for each branch of the diagram. 

The emergence of different processing pathways is a 
consequence of the linear recurrent dynamics: Any activ- 
ity vector x (describing either input s to the network, or 
output r from the network) can be decomposed in terms 
of a sum x = xg +xm- The parting = (x) is a pure baseline 
vector, representing the mean response rate of each neu- 
ron across all stimuli. The remaining part Xm = x — xg 
is a pure modulation vector, with zero baseline. If the 
input is processed by a recurrent network that operates 
linearly on the input according to r = As for some effec- 
tive matrix A (cf. Eq. (3)) it is evident that rg = Asg 
and ~tm = Asm (see Section Theoretical analysis for fur- 
ther explanation). In particular, there is no cross-talk in 



the processing of baseline, sg, and modulation, sm> what- 
soever. The independent, non-interfering processing of 
the baseline and the modulation component of the input 
is exactly corresponding to the two separate pathways 
depicted in Figure 4B. 

For the networks considered in this work, mean and 
variance of all inputs are identical, and all neurons in 
the recurrent network have the same number of exci- 
tatory and inhibitory recurrent afferents. Therefore, the 
entries of the baseline vector sg of the input are all iden- 
tical, and it is mapped to a baseline output firing rate rg, 
which is again a uniform vector. This means that uni- 
form vectors are eigenvectors of the matrices W and A, 
respectively. The eigenvalue belonging to these eigenvec- 
tors is exactly the feedback gain fig of the baseline. In 
networks with dominant inhibition fig is negative. As a 
consequence, the corresponding closed-loop gain yb is a 
positive number, smaller than 1 (Figure 4C). The effec- 
tive enhancement of feature specificity mediated by the 
recurrent operation of the network is the result of differ- 
ent closed-loop gains for modulation, ym, and baseline, 
Yb- For the example network of Figures 1, 2 and 3, this 
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Figure 4 Selective processing of baseline and modulation. (A) General reduced circuit model for the operation of a network on its inputs. 
(B) Reduced circuit model for selective operation of the network on baseline and modulation components of an input vector. (C)Top: Eigenvalue 
distribution of the weight matrix, W, shown for EPSP = 0.1 mV (J = r 5yn eEPSP = 0.136 mV). For normalization, each entry is divided by the reset 
potential, V reS et = 20 mV. The 'exceptional eigenvalue' (green) corresponds to the uniform eigenmode, i.e. the baseline, and the bulk of the 
spectrum (orange) determines the response of the network to perturbations of a uniform input. Bottom, left: Eigenvalue distribution for the matrix 
(1 — W)~ ] , which gives the stationary firing rates. Bottom, right: Sorted magnitudes of eigenvalues of W and (1 — W)~ ] . (D) Baseline and 
modulation components for individual neurons in the network with EPSP = 0.1 mV. Scatter plot (center) shows the modulation vs. baseline 
component of output tuning curves for all neurons of the network. Baseline and modulation are taken as the mean (FO) and the second Fourier 
component (F2) of individual tuning curves, respectively. The markers (cyan crosses) show the center of mass of each cloud. The histograms indicate 
the marginal distributions of baseline (green, top) and modulation (orange, right) components, respectively. 



leads to comparable strengths of baseline and modulation 
in the output tuning curves (Figure 4D), despite having 
much weaker modulation in the input. 

To see how these gains change with the strength of 
recurrent coupling, we fixed the network connectivity and 
only changed the post-synaptic amplitudes in the net- 
work. We then computed the mean baseline and modula- 
tion gains in each network, corresponding to cross marks 
in Figure 4D. Note that the gains are now computed from 
individual tuning curves, ri(0). For each output tuning 
curve, the baseline and modulation component is com- 
puted as the zeroth and the second Fourier component, 
respectively, and then the average values are computed 
across the population (see Methods for details). For the 
cosine tuning we are using here these gains are equal 
to population gains we described before (for rate vectors 
over the population). 

The normalized gains, with respect to corresponding 
gains at zero recurrence, are plotted in Figure 5A. Increas- 
ing the strength of recurrent couplings in the network 
results in a monotonic decrease of yb (green curve), 
since the network is inhibition dominated. The associated 
change of ym (orange curve), however, is non-monotonic: 
It increases until a certain degree of recurrence is reached, 



and then it decreases slowly, while always remaining sig- 
nificantly larger than yb- As a result, the modulation ratio, 
Ym/Yb> of the network exhibits a peak for some degree 
of recurrence (Figure 5B), reflecting optimal performance 
with regard to tuning amplification. 

Stability of the network 

One puzzling observation here is that the behavior of 
networks does not exhibit any dynamic instability as the 
recurrent coupling in the network is increased. This is 
counter-intuitive, as the radius of the bulk spectrum of the 
network, p, increases linearly with the EPSP amplitude in 
the network (see Eq. 29 in Methods). In fact, already at 
EPSP = 0.2 mV this radius is larger than 1, as illustrated 
in Figure 6A, which could result in instability upon stimu- 
lating the network with a modulated input. However, the 
network does not show such instability. 

First, the network dynamics does not change signifi- 
cantly, as demonstrated in Figure 6B and 6C. Moreover, 
the same functional properties are implied from the out- 
put tuning curves of neurons in the network (Figure 7), 
without any sign of unstable operation. Altogether, it 
seems that the networks experience a smooth transition 
as EPSP amplitudes increase, and not an abrupt change of 
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EPSP (mV) 

Figure 5 Network gains for the selective processing of baseline and modulation. (A) Normalized baseline gain, y B norm (green), and modulation 
gain, y^ orr " (orange), for a network with a fixed connectivity, but different strengths of recurrent synaptic couplings. Shaded regions represent 
mean ± std. Lighter colors correspond to lower contrasts. (B) As a result of selective attenuation of the baseline, the normalized modulation-to- 
baseline ratio M orm /yD ,orm ) is generally much larger than I.The degree of recurrence of the network presented in Figures 1, 2, 3 and 4 is marked by 
dashed lines. 



the state. Instabilities are avoided by some dynamic mech- 
anism. We come back to this point in Section Modulation 
gain depends on the operating point of the network. 

A simplified mean-field analysis of the network 

To compute the gains for baseline and modulation of 
inputs, respectively, we employ a simplified mean-field 
approximation, which considers the corresponding aver- 
age gains. For that, we need to compute the mean base- 
line and modulation rate of output tuning curves in the 
network. 

To this end, we first compute the mean firing rate of 
output tuning curves, rg. This is obtained by assuming no 
modulation in the input, i.e. as if all neurons are receiving 
the untuned component of the input tuning curve. This 
is justified by the fact that, in absence of strong non- 
linearities, the two pathways discussed above do not interfere 
with each other. We therefore employ mean field theory 
(see Section Self-consistent firing rates in homogeneous 
random networks) and compute the self-consistent base- 
line firing rates as done previously (Brunei 2000). The 



predicted result, for networks with different couplings and 
for stimuli with different contrasts, is shown in Figure 8A, 
which matches the simulated results quite well. 

Next, we compute the mean modulation compo- 
nent of output tuning curves, rjyi- Here we make an 
approximation: We neglect the modulation of other neu- 
rons in the network and consider only the modulation of 
input to one neuron (see Section Modulation of the firing 
rate and of the membrane potential). This is equivalent to 
assuming 'perfect balance' in terms of modulation, where 
all the modulation components of recurrent inputs from 
the network cancel each other perfectly (/Jm = 0), such 
that only feedforward modulation remains. 

Based on this simplification, the mean modulation of 
the response, tm< is already well predicted (Figure 8B). 
The residual small discrepancy as compared to numerical 
simulations, which is most pronounced for intermedi- 
ate recurrence and high contrast, should be accounted 
for by including network interactions (Section Mixing 
of preferred features in random recurrent networks) 
that amplify the modulation, and spike correlations 
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Figure 6 Stable dynamics of a network with an unstable eingenvalue spectrum. (A) Eigenvalue distribution of the weight matrix, W, shown 
for EPSP = 0.2 mV (J = 0.27 mV). The same normalization as in Figure 4C is employed, i.e. each entry is divided by the reset potential, l/ rese t = 20 mV. 
(B-C) Same plots as in Figure 1 A, B, for EPSP = 0.2 mV. 





Figure 7 Qualitativly similiar functional properties of a network with an unstable spectrum. (A-B) Same distributions of OSI and OSI* as in 
Figure 3A, B, for a network with EPSP = 0.2 mV. (C) Contrast invariance of the OSI measure, as in Figure 3C, for EPSP = 0.2 mV. (D) Distribution of 
baseline (FO) and modulation (F2) component in a network with EPSP = 0.2 mV. Labeling is the same as in Figure 4D. 
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that are ignored in the simplified treatment presented 
here. 

From this result, we can also verify the non-interference 
property with regard to baseline and modulation. The 
separation of pathways proposed in Figure 4B suggests 
that the modulated component of the input vector, sm> 
does not affect the baseline component of output tuning 
curves, rg. This we can verify directly by plotting the mean 
and standard deviation of baseline and modulation com- 
ponents of output tuning curves (Figure 9A). Although 
the variance of the modulation component increases by 
increasing the recurrence in the network, the variance of 
baseline firing rate is very small and does not change with 
recurrence. It only begins to increase when the mean value 
of baseline and modulation components become compa- 
rable in size (about EPSP = 0.1 mV). This is the point at 
which tuning curves experience partial rectification: For 
the tuning curves with modulation component larger than 
the baseline component, some rectification is implied. 
Rectification, in turn, distorts the baseline component of 
tuning curves. 

The non-interference property of baseline and modu- 
lation can also be directly demonstrated from the weight 
matrix. The fact that baseline input does not have any 



component along the modulation vectors became clear 
from the eigenvector of W that corresponds to the excep- 
tional eigenvalue. To show the opposite, namely that an 
input modulation vector induces none, or only negligible, 
baseline in the output response, explicit numerical simu- 
lations of the result of Wsm are performed. The result is 
shown in Figure 9B, which demonstrates that the expected 
value of Wsm (over orientation) is exactly zero, therefore 
not introducing any baseline component, as we discussed 
above. 

Therefore, baseline and modulation are processed sepa- 
rately and independently, with no cross-talk involved, pro- 
vided the network acts linearly on its inputs. In contrast, 
we currently cannot mathematically justify the assump- 
tion of perfect balance. The reason is that the modulation 
vectors, unlike the baseline vector, are not eigenvectors 
of the weight matrix, W. As a result, it is not justi- 
fied to replace W in the product W~tm with a scalar 
value /3m- Treating the problem more rigorously could 
involve expanding the modulation vector in terms of the 
eigenvectors corresponding to the bulk eigenvalues of 
W (Figure 4C), and obtaining the gain accordingly. This 
gain would not, in general, be a single scalar value, nor 
would it be exactly zero, as we have assumed here. We 
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Figure 9 Non-interference of baseline and modulation. (A) Mean and standard deviation of baseline (green) and modulation (orange) 
components for the medium contrast are plotted in logarithmic scales for comparison. As in Figure 8, dots and solid lines indicate simulated and 
predicted values, respectively. Shadings indicate mean ± std. As in Figure 5, mean and standard deviation are evaluated overall neurons in the 
network. (B) The product of the weight matrix W with a baseline, sg, and a modulation input vector, Sm- The entries of the baseline vector are all 
normalized to one, i.e. the input to each neuron is I.The operation of the weight matrix on the baseline vector, Wsg, is plotted in green. The 
operation of the weight matrix on the modulation vector, Wsm, for three different orientations (8 = 30°, 90°, 120°) are plotted in magenta, black and 
cyan, respectively. The corresponding distributions of individual responses are plotted in the magnified histograms on the right. Note that the 
assumption of 'perfect balance' implies a very narrow distribution around zero. The average (over 12 orientations) of the responses, (Wsm)b, are 
plotted in orange. For 12 sample neurons from the network the response vs. orientation of the stimulus are plotted in the inset (center). 



come back to a more precise treatment of the problem in 
Section Linear tuning in recurrent networks. 

Tuning of recurrent inputs 

The assumption of perfect balance of modulation is the 
first-order approximation we make to obtain average gains 
of the network. Here we numerically check how far 
this assumption is from the result of our simulations. 
To answer this question, we investigate tuning of differ- 
ent components of the input to neurons in a network. 
We reconstruct the input from excitatory and inhibitory 
presynaptic sources by replacing each spike with the 
synaptic kernel (alpha function) and adding all the con- 
tributions to obtain a shot-noise signal. We then compute 
the mean value of this signal as the mean presynaptic 
excitation and inhibition, respectively. 

Figure 10A shows these inputs for four sample neu- 
rons from the network. The tuning of excitatory input is 
very weak for all samples; the tuning of inhibition, how- 
ever, is on average stronger. The output tuning of a cell 
results from a combination of contributions from feed- 
forward and from recurrent inputs. The recurrent input 
can therefore change the feedforward tuning: It can either 
retain or change the preferred orientation of the feedfor- 
ward input (upper and lower panels, respectively), and it 
can either amplify or attenuate the tuning strength (upper 
right and lower right, respectively). On average, how- 
ever, the initial shape of the tuning curve is maintained 
(Figure 10B, traces and averages for 12 cells), although 
the tuning of recurrent input leads to a deviation from 
the feedforward tuning, which creates a distribution of 
selectivity. 



Linear tuning in recurrent networks 

The simplified mean-field analysis discussed above 
accounts for the average tuning curve and the mean selec- 
tivity in the network. It does not, however, account for the 
distribution of orientation selectivity across neurons. In 
this section, we resort to a linear analysis of modulation 
processing, in order to provide an approximative analytic 
treatment of this distribution. 

For this linear analysis we need to make two additional 
assumptions. First, we assume that modulations in the 
network can be treated as small perturbations about the 
baseline, and that the dynamics can be linearized about 
this operating point. Note that this assumption implies 
that the contribution of nonlinearities like rectification is 
negligible. Second, we assume that the mixture of tun- 
ings is linear. This assumption is justified as we have 
used cosine tuning (i.e. linear tuning) in the inputs (see 
Methods for details). This allows us to represent each 
tuning curve as a 2D feature vector. Since the operation 
of network on feature vectors is linear, the mixture of 
tunings is now reduced to the vectorial summation of 
corresponding tuning vectors (see below, and Methods). 

We therefore start first by linearizing the dynamics 
about the operating point, i.e. the baseline. After com- 
puting the baseline firing rate, rg, as described before, we 
compute the mean and standard deviation of the input to 
each neuron in the baseline state (/zg and a%, respectively; 
Eq. (32) in Methods) as a function of baseline input sg and 
baseline firing rate rg 
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Figure 10 Tuning of neuronal inputs. (A) Tuning of the mean input from excitatory (red) and inhibitory (blue) neuron populations, and of the 
external input (green). Shown are four sample neurons, for stimulations at the medium contrast. EPSP = 0.2 mV. The total recurrent input 
(excitatory + inhibitory), and the net modulation of the input (external + total recurrent) are also plotted (in black and orange, respectively). The 
input is computed by replacing each presynaptic spike by an alpha kernel and computing the mean amplitude of the shot-noise signal (in mV/ms). 
1 mV/ms corresponds to a mean membrane potential at the spike threshold of the neurons in our simulations. (B) Same for twelve sample neurons, 
along with their mean values (thicker lines). 



We then compute the linear gains by perturbing the input 
with a small Ss. We dismiss all contributions of order two 
or higher in the perturbation parameter and write 

Sr = UsSs, (4) 

where Sr is the change in the output firing rate result- 
ing from a perturbation of the input rate Ss. Here, we 
compute these gains numerically by solving the perturbed 
mean field equations about the basline (see Methods, and 
Eq. 12). Figure 11A shows the input-output relationship 
for the network with EPSP = 0.2 mV at all contrasts. 

After computing the linear gains, we can rewrite the 
linear rate equations of the network for the modulation 
pathway. The modulation in the firing rate of each neu- 
ron is caused by a linear mixture of external input and the 
input it receives from the recurrent network 

r M = S(Wr M + J s sm)- (5) 

Both of these terms are now weighted with the extra 
factor, f, from linearization, which is the effective gain 
of modulation at this operating point. If we rewrite the 
equation as 

00 

m = (i - 1 wr^hsM = a w) k s M , (6) 

k=0 

it can further be approximated by 

r M » SJs(sm + SWsm). (7) 



Here we are neglecting the contribution of higher- 
order recurrent inputs ((^W) 2 sm> (? W) 3 sm, ■ ■ ■) in the 
processing. 

Next step is to interpret the above equation for tun- 
ing curves. Since we have assumed cosine tuning for the 
input, we can represent each input tuning curve by a vec- 
tor, Si. Likewise, we represent output tuning curves by 
vectors Rt. We refer to these vectors as the Tuning Vec- 
tors. For a 2D feature like orientation selectivity we obtain 
2D Tuning Vectors. The angle of this vector represents 
the input preferred orientation 0*, and the length of it 
is a measure of orientation selectivity (it is indeed equal 
to OSI before normalization). For notational convenience, 
we represent these 2D vectors by complex numbers. We 
identify real parts with x-directions, and imaginary parts 
with y-directions. Any 2D vector then corresponds to a 
complex number in a one to one fashion. 

The mixture of tuning curves in cosine tuning is now 
simplified to the summation of the corresponding Tun- 
ing Vectors in the complex plane (see Methods for further 
details). We can therefore rewrite the linear rate equation 
Eq. (7) for Tuning Vectors 

R = SJ s (S + t.WS). (8) 

Here R and S are vectors with complex elements, rep- 
resenting output and input tuning of all neurons in the 
network, respectively. This is now an equation which 
expresses the output Tuning Vectors in terms of a linear 
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Figure 1 1 Linear tuning in recurrent networks. (A) Linearized gains for single neurons embedded in the network. The extra firing rate, &r, of a 
neuron produced in response to a small perturbation, J 5 Ss, in the input intensity, plotted for different baseline inputs corresponding to different 
contrasts. The response is computed by numerically perturbing the mean-field equations (see Methods). The linear gain, J = Sr/(J s Ss), is then 
computed by linear regression of data points. For this example with EPSP = 0.2 mV, the value f = 0.026 is obtained, which is also used for the next 
panels. (B) For the sample neuron in Figure 2A, all presynaptic Tuning Vectors are extracted (Vexp(/2#*), weighted by the linear gain (J), and 
vectorially added together, reflecting linear integration in neurons. Although each presynaptic vector makes only a small contribution, the resulting 
random sum can lead to a large resultant Tuning Vector. These are generally larger for presynaptic inhibition (Presyn. Inh.) compared to presynaptic 
excitation (Presyn. Exc). Note different scales of axes. (C) Left panel: The resultant vectors for recurrent excitation (Rec. Exc), recurrent inhibition 
(Rec. Inh.), total recurrent (Rec. Tot. = Rec. Exc + Rec Inh.), feedforward input (Input), and the total input (Tot. = Input + Rec Tot.) are plotted. All 
normalized input Tuning Vectors have the same length of one, denoted by the green circle. Right panel: Total recurrent Tuning Vectors (Rec. Tot.) for 
all neurons in the network are compared with the normalized length of their input Tuning Vectors (green circle). D) Distribution of the length of all 
Tuning Vectors for all the neurons in the network. Dashed lines show the predicted distributions of the linear analysis in each case (see text for details). 



mixture of input Tuning Vectors. Similar to Eq. (7), we 
are neglecting higher-order terms ((f W) 2 S, (f W) S, . . .) 
here. 

An individual output Tuning Vector, 7?/, is then given 
by Ri = t;f s (Si + t, J^j WySj): The output tuning of each 
neuron is a mixture of its input tuning and weighted vec- 
tors of all presynaptic sources. For the specific example of 
the neuron in Figure 2A, all the contributions of presy- 
naptic sources are shown in Figure 11B, for excitatory 



and inhibitory populations separately. Each small jump 
in the space represents the contribution of a presynaptic 
Tuning Vector, the size of which is f/ or — £gj for excita- 
tory and inhibitory presynaptic sources, respectively. Si is 
normalized to 1. 

Although each presynaptic contribution is much smaller 
than the input (of order £/ = O(10 -3 ) in this case), 
the resultant vector (dashed lines) can be large. In par- 
ticular, the resultant vector of inhibition is comparable to 
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the length of the input vector for this specific example 
(Figure 11C, left). Since the angle of the vector is close 
to the input angle, this leads to an amplification of the 
resultant tuning, although it typically also changes the 
preferred orientation (Figure 11C, left). This explains why 
the OSI of this neuron (0.65) is larger than the mean OSI 
of the network (0.42, see Figures 2 and 3). 

Not all the vectors resulting from recurrent contri- 
butions, however, are large, nor do all have a similar 
preferred orientation as the input tuning. Indeed, as the 
connectivity is random and the initial preferred orienta- 
tions are assigned randomly to each neuron, the preferred 
orientation of the resultant vectors are also random. This 
is shown in Figure 11C, right panel, where all the recur- 
rent Tuning Vectors are explicitly computed (by reading 
the input Tuning Vectors and the connectivity), and plot- 
ted in the complex plane. The distribution of the vectors in 
this plane is a 2D Gaussian distribution, according to the 
Central Limit Theorem. As a result, most of the Tuning 
Vectors have small magnitude, and only a few of greater 
magnitude contribute to the tail of the distribution. 

The distribution of this length is plotted in Figure 11D, 
for different subpopulations of neurons. The peak of the 
distribution for excitatory neurons is at smaller values 
than for inhibitory neurons, and the overall length of 
recurrent tuning is mainly determined by inhibition in the 
network. Knowing the standard deviation of distributions 
of Tuning Vectors, one obtains the distribution of vector 
lengths according to ^e~ x /(2fr ' (see Methods, Eq. (62)). 
In this example, the standard deviations are er e xc = 0.11, 
cr in h = 0.44 and cr t ot = 0.45, for excitation, for inhi- 
bition, and for all recurrent neurons, respectively. The 
length of tuning vectors predicted by this result is plot- 
ted in Figure 11D (dashed lines). The length of overall 
tuning, i.e. recurrent Tuning Vectors vectorially combined 
with the input Tuning Vector (green), can also be com- 
puted (see Methods, Eq. (63)). Normalizing the input 
Tuning Vectors to length 1, the distribution amounts to 
^ e -(* 2 +i)/(2 ( r 2 ) /o( ^ )| p i otted in the same figure (dashed 
purple line). 

From Figure 11D one can compare the strength of 
the input tuning (green line) with the tuning gener- 
ated within the random network (black distribution). The 
mean length of the recurrent tuning is smaller than the 
feedforward, and only few neurons show comparable tun- 
ing strength. The distribution of the combined tuning 
strength (purple line), which is a mixture of feedforward 
and recurrent components, has now a broad distribu- 
tion, where many neurons show less tuning than their 
input (less than 1, attenuated), and many more have 
enhanced selectivity (greater than 1, amplified). In gen- 
eral, amplification happens when the randomly generated 
Tuning Vector within the network is roughly aligned with 
the initial Tuning Vector, and attenuation happens for 



recurrent Tuning Vectors in the opposite directions' 5 . A 
random recurrent network, thus, in itself is capable of 
attenuating and amplifying orientation selectivity. This is 
a mechanism in addition to the selective gains of base- 
line and modulation described before. This mechanism, 
however, comes at the expense of shifting the tuning 
curves of neurons from their initial, feedforward preferred 
orientations. 

Regimes of orientation selectivity 

As Figure 12A shows, in a weakly coupled network the 
PO of neurons are hardly changed with respect to their 
input PO. This is a regime where feedforward projections 
are dominant with respect to functional properties of 
neurons. Under these conditions, the recurrent network 
cannot compensate the increase in the baseline, and both 
baseline and modulation are scaled identically (Figure 13). 
As a result, neuronal tuning curves tend to simply reflect 
the tuning of the input (Figure 13), and tuning curves 
reduce their selectivity for high-contrast inputs, because 
the feedback compensation is weak. Although the average 
orientation selectivity index (OSI) of neurons in the net- 
work could be high for the lowest contrast in a weakly con- 
nected network, this selectivity is lost when the contrast is 
increasing (Figure 14A). 

As the strength of recurrent couplings increases, the 
contribution of the network becomes more effective to 
attenuate the baseline and selectively enhance the modu- 
lation (Figure 13 and Figure 5). This leads to more stable 
OSI distributions across different contrasts (Figure 14B 
and 14C), and hence makes the selectivity more robust. 
Moreover, as a consequence of stronger recurrence in the 
network, output POs deviate more from their initial PO 
(Figure 12B and 12C), since the strength of recurrent con- 
tributions (recurrent Tuning Vectors) has now increased. 
This is summarized for all networks in Figure 12D by a 
scatter degree index (SDI), which quantifies the degree of 
PO deviation in each network. 

Notably, SDI does not linearly increase with recurrent 
connection strength. Rather, it saturates for rather weak 
connections, and then reaches an asymptotic value. The 
reason for this behavior is that the contribution of recur- 
rent Tuning Vectors in the final tuning depends on / and 
the linear gain, f , as we described in the previous section 
(see Eq. (8)). Although / is monotonically increasing in 
Figure 12D, the effective strength of recurrent Tuning 
Vectors depends on the product /£ . It appears, therefore, 
that the linear gains are not increasing as / increases, or 
they are even decreasing. 

This trend is further supported by shape and size of 
the tuning curves (Figure 13). For networks with strong 
recurrent coupling, the maximum firing rate of tun- 
ing curves decreases and the modulation component 
becomes smaller. Since the linear gains determine the 
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Figure 12 Input vs. output preferred orientation of neurons in the network. (A-C) Output PO vs. Input PO for three values of recurrence, at the 
medium contrast. Increasing the recurrence leads to more scatter about the diagonal. For illustration purposes, the excitatory (red) and inhibitory 
(blue) neurons have been plotted only above or below the diagonal, respectively. (D)To quantify the amount of PO change when going from input 
to output, the Scatter Degree Index (SDI) is plotted as the angular deviation of APO = OutputPO — InputPO (see Methods). The maximum value of 
this index is 40.5°, which corresponds to a uniform distribution of APO. Darker colors show higher contrasts, respectively. 



embedded gain of neurons in the network in response 
to modulations, this trend also suggests that these gains 
are decreasing in the highly recurrent regime. This was 
indeed visible in the behavior of gains and firing rates for 
modulations, shown in Figures 5A and 8A, respectively. 

Although increasing the recurrence stabilizes the OSI, 
it makes the neurons of the network less feature selec- 
tive, if the recurrence is too large (Figure 14A-D). There 
is, therefore, a trade-off between orientation selectiv- 
ity and contrast invariance in the networks. Increasing 
the recurrence makes the negative feedback in the base- 
line pathway stronger, making the divisive suppression 
of the baseline - and hence the contrast invariance - 
more effective. This comes, however, at the expense of a 
decrease in the gain in the modulation pathway, which 
makes the responses weaker and less selective. We have 
quantified this trade-off by dividing the mean OSI of indi- 
vidual tuning curves in a network by its standard devia- 
tion across different contrasts. The intermediate recurrent 
regime shows optimal behavior with large and stable OSI 
(Figure 14D), and it more or less coincides with the region 
of optimal tuning amplification (Figure 5). 

Tuning and invariance of membrane potentials 

It is also informative to look at the membrane potential 
dynamics of neurons in the network. Figure 15A shows 



the membrane potential of a sample neuron in response to 
a stimulus of its preferred orientation, and the orthogonal 
one. The mean membrane potential remains, on average, 
below threshold, as it is reflected in the distribution of 
membrane potential (Figure 15A, right panel), and only 
fluctuations in the input leads to spiking activity. For the 
orthogonal orientation, this activity is very sparse; for the 
preferred orientation it leads to a reasonable firing rate. 
If we plot the mean membrane potential for each orien- 
tation, it shows the same tuning as of the spiking activity 
(Figure 15B). 

Plotting the free membrane potential for this neuron 
(Figure 16A, left), and more samples from the same 
network (Figure 16A, center), verifies this observation. 
The free membrane potential (Figure 16A, right) remains 
below threshold and has a similar tuning. Indeed, the tun- 
ing is slightly enhanced, since the negative contribution 
of the reset voltage is now corrected for. This behav- 
ior is consistent across different contrasts. Increasing the 
contrast shifts the mean tuning curve as a whole down 
to more negative values, as a result of more negative 
feedback recruitment in the network. This, in turn, com- 
pensates for the increase in baseline, and suppresses the 
response to non-preferred orientations. The tuned part, 
however, is still capable of generating spikes, which leads 
to the tuning of spiking activity. 
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Figure 13 Neuronal tuning curves in weakly (left) and strongly (right) recurrent networks. (A) Tuning curves of a sample neuron (same as in 
Figure 2A), for different degrees of recurrence in the network, as indicated by different EPSP amplitudes. (B) More (125) sample tuning curves from 
the network, aligned to their Input PO. Red and blue curves show excitatory and inhibitory output tuning curves for the medium contrast, 
respectively. Shown in green is the input tuning curve, normalized to the average (over the population) of the mean (overall orientations) of all 
tuning curves in the network. (C) Mean and standard deviation (across neurons) of all aligned tuning curves, for networks with different degrees of 
recurrence. Lighter shadings denote lower contrasts, respectively. 



If the recurrent compensation was not effective, a dif- 
ferent behavior would emerge. In fact, if the recurrent 
coupling is very weak (Figure 16B, left), the free mem- 
brane potential is above threshold for almost all orien- 
tations. In such a case, the response to all orientations 
is in the mean-driven regime, which yields significant 
firing rates for the preferred as well as orthogonal orien- 
tations. As a result, the so-called iceberg effect broadens 
the tuning curves significantly upon increasing the con- 
trast. Increasing the recurrent coupling shifts the mean 
membrane potential down and the network operates in 
the fluctuation driven regime; this makes the tuning of 
membrane potential and the resulting spiking activity 
more robust and contrast invariant (Figure 16B, center 
and right panels). Indeed, for the intermediate recurrent 
regime (Figure 16B, center) the tuning is perfectly contrast 
invariant. 

Spiking activity in inhibition-dominated networks 

The recurrent excitation in inhibition-dominated net- 
works of the sort we are studying here is over- 
compensated by the surplus recurrent inhibition. Some 
residual inhibition remains as the net effect of recurrent 
interactions. If the recurrent coupling is strong enough, 
this net inhibition determines the effective threshold of 



neurons in the network. Therefore, it is not the thresh- 
old mechanism of neurons which cuts off the responses 
at non-preferred orientations. Balance of excitation and 
inhibition within the network, in contrast, governs the 
generation of spikes and, hence, attenuation and ampli- 
fication of the baseline and modulation components, 
respectively. 

A sample of this temporal balance is shown in 
Figure 17A for the net excitatory and inhibitory inputs 
from the network to a neuron. Although the net excitation 
is on average above the firing threshold, the net inhibition 
is twice as large on average, as a result of the parameter 
configuration used. Moreover, occasional deflections in 
excitation and inhibition follow each other on a fine time 
scale. The net recurrent input to the neuron is on average 
negative. Occasional disbalance, however, provides a net 
excitatory drive for brief periods of time. 

Considering spike-triggered averages of the net excita- 
tory and inhibitory input (Figure 17B) reveals that spike 
emission in the network is mainly governed by recurrent 
inhibition, rather than recurrent excitation. Therefore, in 
the inhibition-dominated regime, fluctuations of inhibi- 
tion is the main determining factor for spiking activity, 
in agreement with the results of experimental studies 
(Rudolph et al. 2007). 
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Figure 14 Distribution and contrast invariance of selectivity. (A-C) Distribution of orientation selectivity in three different regimes of 
recurrence. Lighter colors code for lower contrasts. The mean OSI for each distribution is indicated in the corresponding brightness, respectively. 
OSI* is computed from the cosine fit, as in the main panel of Figure 3B. (D) The mean OSI of all neurons in the network is shown for different levels 
of recurrence, at three different contrasts. Lower contrasts are plotted in lighter colors. Increasing the recurrence makes the OSI less susceptible to 
changes in contrast. For high recurrences, this invariance comes at the expense of a decreased selectivity, as the mean OSI in the network decreases. 
This trade-off between selectivity and invariance is quantified (in brown) by the average (over contrasts) of mean OSI (across neurons) divided by 
the standard deviation (over contrasts) of the average OSI (across neurons). 
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Figure 15 Tuning of the membrane potential. (A) Membrane potential of a sample excitatory neuron from the same network with 
EPSP = 0.2 mV, at the medium contrast. Traces of the membrane potential are plotted for 1 s of response to the preferred (red) and its orthogonal 
(cyan) stimulus orientation. The histograms of the membrane potential for 6 s of stimulation are shown on the right. (B) Traces of the membrane 
potential along with the elicited spikes for 12 orientations, for 1 s of recording. The tuning curves of the mean membrane potential (red) and the 
corresponding firing rate (black) is computed from 6 s of stimulation in the right panel. 
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Figure 16 Membrane potential at different contrasts. (A) Left: Tuning of the membrane potential of the same neuron as in Figure 15 for 
different contrasts. Center: Average tuning curve of 24 (12 excitatory, 12 inhibitory) randomly sampled neurons, ranging over all POs. The tuning 
curves are aligned at a PO of 90°. Error bars indicate mean ± std over sampled neurons. To improve the display, they are plotted only for every third 
data point. The mean membrane potential (over all neurons and all orientations) is indicated by solid, horizontal lines in each case. The shading 
represents the standard deviation (over neurons) of the mean (over orientations) of the average membrane potential Wb). Right: Same as the panel 
in the center, for the free membrane potential, Vf m = V + rrV !eset . (B) The same mean ± std (over sampled neurons) of tuning curves for the 
distance to threshold of free membrane potentials (Vf m — l/ t h), for different regimes of recurrence. For the lowest contrast, most of the error bars are 
smaller than the size of markers and hence not visible. The mean membrane potential and its standard deviation over the (sampled) population is 
shown, as before, by horizontal lines and shadings, respectively. 



Modulation gain depends on the operating point of the 
network 

The strength of net inhibitory feedback also affects the 
selective gains. Whereas the mean inhibitory recurrent 
input sets the divisive gain in the baseline pathway 
(Figure 4B), it affects the modulation gain by determining 
the mean distance of the membrane potential to thresh- 
old. As suggested by Figure 16B, increasing the level of 
recurrence induces larger average distances of the mean 
membrane potential from threshold. This is shown for all 
networks in Figure 18A. 

Using the simplified mean-field analysis provided in 
Section A simplified mean-field analysis of the network, 
we can predict the mean membrane potential of a net- 
work. Knowing the baseline firing rate of the network, r B , 
the mean baseline membrane potential of the network, vg, 
is obtained (see Eq. (25) in Methods) as 

vb = r [-Vthrs + fiB] > (9) 



where \x B = JNe(f — g(l — f))rg + J s sg (see Methods, 
Eq. (32)). Note that, as the input is shunted during the 
refractory period, the shunted feedforward and recur- 
rent input should be subtracted from the total input 
in this expression. The shunted input can be computed 
as r B t le [fiB (Kuhn et al. 2003; Murthy and Fetz 1994). 
The corrected membrane potential, after considering the 
effect of refractory period, is then given as 

v B = r [-Vthrs + (1 - r B t re i)lJ-B] • (10) 

The result of this prediction for different contrasts is 
plotted in solid lines in Figure 18A. 

Note that this prediction is obtained under the assump- 
tion of a Gaussian distribution of input to all neurons. The 
prediction, therefore, fails to be exact if this assumption is 
violated. The distribution is, in fact, skewed as a result of 
correlations in the network (Kuhn et al. 2003). The devi- 
ation from a Gaussian distribution increases for higher 
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Figure 17 Spiking activity in the inhibition dominated regime. (A) Statistical balance of overall excitatory (red) and inhibitory (blue) input from 
the recurrent network. The dominant inhibition keeps the net recurrent input (gray) negative on average, and only occasional transients lead for a 
moment to a net excitatory drive from the network. EPSP = 0.2 mV. (B) Spike triggered averages (STA) of excitatory (Exc.) and inhibitory (Inh.) 
recurrent inputs are plotted from spikes of 1 2 randomly sampled neurons in response to one stimulation (6 s). The total recurrent input (Tot.) is 
plotted in black. The membrane potential is normalized by V t h (Norm. Vm.). 
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Figure 18 Mean distance to threshold sets the operating point of the network. (A) The mean membrane potential for networks with different 
EPSP sizes (circles), at different contrasts, along with the predicted values (solid lines). Mean membrane potential is computed as the average (over 
orientations) of mean tuning curves of 12 sample neurons (the same as in Figure 16A). Lighter colors belong to lower contrasts, respectively. 
(B) Input modulation normalized by the distance to threshold, Vth — Vm, (solid lines) compared to the output modulation (orange circles; same as 
Figure 8B). Input modulation is given as the input modulation rate times its efficacy (J s mss). 
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recurrences, explaining the discrepancy of our predictions 
for highly recurrent regimes. 

The mean membrane potential is crucial in determin- 
ing the overall gain of the linearized dynamics. A very 
depolarized membrane potential reduces the chance of an 
input perturbation to reach the firing threshold and to 
elicit a spike. Therefore, it affects the gain of modulation, 
as shown in Figure 18B: The mean output modulation of 
tuning curves, ym> is indeed inversely related to the mean 
distance to threshold. This suggests that the embedded 
gain of neurons in the network in response to modulation 
is also inversely proportional to the distance to thresh- 
old. Although the mean modulation in the input is below 
the spike threshold of a single neuron, fluctuations are 
nevertheless capable of eliciting reasonable firing rates. 
The resultant linearization of the /-/ curve, as shown 
in Figure 11A, is a result of the noise, er#, generated 
within the network due to the balance of excitation and 
inhibition. 

If we now compute these gains for networks with dif- 
ferent degrees of recurrence, the linearized gains match 
well with the mean modulation gain in the network 
(Figure 19, for the middle contrast). This suggests that the 
mean amplification of modulation in the network is fully 
accounted for by the linear gain, which in turn depends on 
the operating point of the network as defined by the mean, 
[AB> and the standard deviation, erg, of the input to neurons 
from the network in the baseline state. The linear gains are 
obtained by perturbing r = f(i±, er) at this operating point 



df(u,a) df(u,a) 
8r= J — - g/z+ J ^ ho, (11) 
dfl da 

and using / s f = Sr/Ss. Here, we determine this quantity 
numerically as 

SsM = ffais + 8s), a (s + 8s)) -f(fi(s),a(s)) (12) 

for Ss = 100 spikes/s (Figure 19). 

In summary, the inhibitory feedback in a recurrent net- 
work contributes to orientation selectivity in crucial ways. 
First, it provides a negative feedback which offsets the 
baseline component of the input tuning curves and leads 
to divisive attenuation of the common mode (selective 
attenuation in the baseline pathway). Second, it sets the 
operating point of the network and determines the lin- 
earized, embedded gain, which in turn determines the 
modulation gain (selective amplification in the modula- 
tion pathway). Moreover, the feature selectivity generated 
by the recurrent network as a result of summing many 
inputs of random selectivity leads to either amplifica- 
tion, or attenuation, of the feedforward tuning (random 
summation of recurrent tuning). Since the contribution 
of each presynaptic modulation vector must be weighted 
according to the linearized gain, the bulk spectrum of the 
connectivity matrix W must also to be weighted accord- 
ingly. The spectrum of the network with EPSP = 0.2 mV 
shown in Figure 6A, for instance, was obtained by normal- 
izing / by V t h- The linear gain suggests now that / should 
be multiplied by f = 0.026, which is by a factor 2 smaller 
than 1/Vth = 0.05. This leads to a decrease in the radius 



0.08 
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Figure 19 Linear gains determine the gain and stability of networks in response to modulation, or the medium contrast, linear gain (J, black 
line) is computed at each EPSP and is compared with the corresponding modulation gain, = Output modulation/Input modulation (orange), f 
is computed as f = Sr/(J 5 Ss), for a small perturbation of the input, Ss = 100 spikes/s. Inset: The radius, p, of bulk spectrum of W, normalized by the 
linear gain ({ ) at each EPSP. Instead of dividing J by V± (as in Figure 4C), J is now multiplied by f . As a result, the normalized radius is now obtained 

asp norm _ ? jJfaQ _ e ) [f + g2(1 _f)]. 
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of the bulk from p = J/V^Ne(l - 7) [f +g 2 (l -/)] « 

1.68 to p = SI^NeO- - 7) [f + g 2 (l -/)] « 0.87. This 
implies that none of the modulation modes corresponding 
to the bulk are actually unstable, at this operating point 
of the network. Indeed, if we plot the normalized radius, 
p norm , for all the networks at different contrasts, the p norm 
never exceeds one (Figure 19, inset). This means that, 
although the coupling strength is monotonically increas- 
ing, the network dynamics stabilizes the spectrum in 
inhibition-dominated networks (see Pernice et al. (2012) 
for related observations). 

Discussion 

Using large-scale simulations and associated mean-field 
analysis of networks of spiking neurons, we have demon- 
strated how highly selective neuronal responses can be 
obtained in random networks without any spatial or 
feature specific structure. Our mathematical analysis pin- 
points the mechanisms responsible for selective attenu- 
ation of the common mode and selective amplification 
of modulation, and predicts some essential properties of 
these networks. 

A generic model of local circuitry 

Here we discussed the specific case of orientation selec- 
tivity in the early visual system, as we were able to link 
our findings to an ample body of experimental litera- 
ture. However, our model could potentially be of a much 
broader scope. It proposes a general mechanism for the 
emergence of strong feature selectivity, which could actu- 
ally be at work in other sensory modalities as well. Our 
network model can thus be conceived as a generic model 
for the local cortical circuitry, which enhances feature 
selectivity and ensures contrast invariance of process- 
ing, without resorting to feature specific structure or 
experience-dependent fine tuning. 

Our analysis suggests that a randomly connected net- 
work with dominant inhibition is already capable of selec- 
tively removing the uninformative common mode of a 
stimulus that is represented by the network in a dis- 
tributed fashion, while preserving the informative mod- 
ulations in the response pattern induced by stimulation. 
This way, the tuned part gains salience, and the signal- 
to-noise ratio improves. The network also amplifies the 
tuned component (signal) by two mechanisms: First, by 
modulating the embedded gain through adjusting the 
operating point of the network, and second, by recurrently 
mixing presynaptic selectivities and thereby amplifying 
weakly tuned inputs in some neurons. 

Regimes of orientation selectivity 

The same mechanisms could also lead to attenuation of 
the signal in the network. First, increasing the recurrent 



coupling in the network increases the mean distance of the 
membrane potential to the firing threshold, which in turn 
decreases the modulation gain in the network. Second, the 
recurrent mixing of weak tunings in the network gener- 
ates a distribution of selectivity, with attenuation in many 
neurons. As a result, this mechanism cannot increase the 
selectivity of a certain fraction of the neuronal population 
beyond the input selectivity, unless this is compromised 
by a decrease in the selectivity of another fraction of the 
population. 

In addition, there is a trade-off between selectivity and 
contrast invariance of the neuronal responses. Increasing 
the degree of recurrence in the network makes the selec- 
tivity more invariant and the distribution of it more robust 
with regard to variations of contrast, but it decreases the 
overall gain of modulation. As a result, there is a region 
of intermediate recurrence in our networks, where tuning 
amplification is most pronounced, and orientation selec- 
tivity has the largest value with the lowest sensitivity to 
contrast. 

Our computational study, therefore, suggests an inter- 
mediate regime of recurrence as the optimal regime of 
feature selectivity for early sensory processing. This is the 
state of the network which exhibits the stimulus driven 
properties of neurons observed in experiments (Hofer 
et al. 2011), while preserving other important features like 
strong and contrast invariant orientation tuning. In this 
regime, the feature selectivity of neurons would exhibit 
the least deviation from their input selectivities, essen- 
tially reflecting the tuning of the feedforward input. The 
role of the recurrent network at this stage would then be 
to enhance this selectivity, by performing operations like 
increasing the signal-to-noise ratio and contrast invariant 
gain control. 

This scenario might best explain the state of the input 
layer L4, in which orientation selectivity first emerges in 
the cortex. The same is not necessarily true for more 
recurrent layers like L2/3 or L5, which are involved in 
later stages of sensory processing like learning, association 
and motor control. It is, therefore, plausible that differ- 
ent regimes of recurrence exist in different layers, which 
may be suited to perform different types of processing. 
One measure for the degree of recurrence in a network 
is the tuning of the total recurrent input. As the recur- 
rent coupling increases in a network, the tuning of the 
recurrent input generated within the network increases as 
well, and the assumption of untuned total input becomes 
a questionable approximation. 

There are indeed contradictory results reported in 
experimental studies on the tuning of input in rodents: 
Both untuned inhibition (Li et al. 2012a,b; Liu et al. 2011) 
and co-tuning of excitation and inhibition (Tan et al. 201 1) 
have been reported by different laboratories. Our results 
show that even in absence of significant tuning of the total 
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input received from the recurrent network, another mech- 
anism of selective attenuation and amplification can lead 
to strong selectivity and contrast invariance. This, how- 
ever, does not exclude a random summation of selectivity 
within the recurrent network as a contributing mecha- 
nism. Indeed, in the first example network we investigated 
here, both mechanisms were at work. It is, therefore, 
plausible that both tuned and untuned components exist 
to some degree in such networks, but the exact mixing 
depends on the operating point and, specifically, on the 
degree of recurrence. This would therefore suggest that in 
more recurrent layers like L2/3 more neurons with strong 
input tuning should be recorded, while in the input layer 
L4 untuned inputs preponderate 0 . 

It should be noted, however, that our discussion here 
applies to recurrent excitation and inhibition. Tuned exci- 
tation and inhibition, when measured in terms of the 
total excitatory and inhibitory conductances in intracellu- 
lar recordings, are the total excitatory and inhibitory input 
that a neuron observes. It is therefore possible that the 
tuning of the feedforward input is dominant in the tuning. 
Likewise, feedforward inhibition, mediated by disynap- 
tic inhibition, can have the same tuning as the feedfor- 
ward excitatory input, as the former is mediating it. The 
mechanisms discussed here, however, apply to recurrent 
excitation and inhibition, since they are a consequence of 
the dynamics of a network of synaptically connected neu- 
rons and, in particular, recruitment of feedback inhibition 
within the network. 

Recurrent vs. feedforward inhibition 

Recurrent inhibition in our networks selectively feeds the 
mean signal back and subtracts it from the tuning curves. 
The overall effect of this subtraction results in a divi- 
sive attenuation of the baseline. This untuned suppression 
has been experimentally demonstrated to play a crucial 
role in increasing the selectivity (Shapley et al. 2003; Xing 
et al. 2011). More specifically, it has been recently shown 
that in L4 of mouse visual cortex, it underlies sharpen- 
ing of orientation selectivity (Li et al. 2012a,b). This is 
also consistent with the results of a recent study on the 
role of somatostatin expressing, SOM + , GABAergic neu- 
rons in orientation selectivity (Wilson et al. 2012) (but 
see Lee et al. (2012b)). The subtraction of the baseline, 
attributed to this specific subtype of inhibition by Wilson 
et al. (2012), effectively leads to the selective attenua- 
tion/division of the baseline, as described in the present 
article. Reduction of this inhibition would therefore lead 
to a constant increase in the baseline activity of the tuning 
curves, which has indeed been recently reported in Dlxl 
(-/-) mice with selective reduction of activity in dendrite 
targeting inhibitory interneurons (Mao et al. 2012). 

This mechanism is in contrast to the role of feedfor- 
ward inhibition of fast spiking interneurons, parvalbumin 



expressing, PV + , GABAergic neurons. As opposed to 
SOM + neurons, which are more recurrent, these neurons 
are mainly driven by feedforward input (Adesnik et al. 
2012). Unlike SOM + neurons, which are involved in den- 
dritic computation and controlling the input, PV + neu- 
rons are better suited for controlling the output, as they 
innervate the peri-somatic regions (Di Cristo et al. 2004; 
Ma et al. 2010). Consistently, they are also more effective 
during the transient responses, as reflected in their activa- 
tion pattern (Tan et al. 2008). In contrast, SOM + neurons 
are better suited for sustained activity. 

These properties might then suggest that feedforward 
inhibition is primarily involved in gain control (Atallah 
et al. 2012), which uniformly rescales all components of 
the tuning curves. The attenuation of the baseline, there- 
fore, comes at the expense of attenuating the modulation. 
It cannot be selective to the baseline, in contrast to the 
recurrent mechanism we have modeled here. Note that, 
for simplicity, we have not considered feedforward inhibi- 
tion in our model. However, feedforward inhibition could 
easily be added to the model by mediating the same exci- 
tatory input to each neuron by an inhibitory neuron. 
Doing so would effectively lead to a change in the over- 
all feedforward gain, provided that inhibition is not strong 
enough to result in rectification (compare with Lee et al. 
(2012b)). 

If PV + neurons are strongly driven by feedforward 
input, and if the feedforward input is only slightly tuned, 
as we assumed here, the responses of PV + neurons should 
be only weakly modulated. In contrast, as the results of our 
simulation showed here, inhibitory neurons involved in 
recurrent computations can be highly selective, although 
they receive weakly modulated inputs. In agreement with 
this, SOM+ inhibitory neurons involved in recurrent inhi- 
bition have been reported to be as selective as excitatory 
neurons, in contrast to PV + neurons with broader selec- 
tivity (Ma et al. 2010). However, it should be possible 
to make PV + interneurons more selective, by providing 
more recurrent inhibition to them. In fact, it has recently 
been reported that activating SOM+ inhibitory neurons 
can unmask and enhance the selectivity of PV + cells by 
suppressing untuned input (Cottam et al. 2013). 

Note that, as contrast invariance depends on the selec- 
tive attenuation of the baseline, it should be the result of a 
recurrent mechanism. We therefore suggest that the con- 
stant increase of untuned inhibition that neurons receive 
upon increasing the contrast (Li et al. 2012b) should be a 
result of the recurrent, and not of the feedforward, inhi- 
bition. This may explain why dark reared mice in this 
experiment, which lacked a broadening of PV+ responses, 
still show an aggregate untuned input from the network 
and hence highly selective responses (Li et al. 2012b). 
Although individual inhibitory neurons were on average 
highly selective in our networks, the emergent result of the 
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interaction of excitation and inhibition lead to an effec- 
tive untuned inhibition, which increases proportionally 
with contrast. This is again consistent with the results of 
Li et al. (2012b), who demonstrated that "blocking the 
broadening of output responses of individual inhibitory 
neurons does not block the broadening of the aggregate 
inhibitory input to excitatory neurons". It is also consis- 
tent with the results of a previous report, demonstrating 
that "broad inhibitory tuning" of fast spiking cells is "not 
required for developmental sharpening of excitatory tun- 
ing" (Kuhlman et al. 2011). Based on these results, we 
therefore hypothesize that untuned inhibition might be an 
emergent property of an inhibition dominated network, 
and not a feedforward consequence of broadly tuned fast 
spiking neurons. 

Comparison with other models 

Most existing recurrent theories of orientation selectivity 
consider the case of species like carnivores and primates, 
with a clustered organization of selectivity in orienta- 
tion maps (Ben-Yishai et al. 1995; Somers et al. 1995). 
Consistent with the proximity of neurons with similar 
selectivity, these theories assume a feature specific con- 
nectivity of neurons. The Mexican hat profile of con- 
nectivity which they assume leads to a more broadly 
tuned inhibition, which suppresses the mean, and to a 
sharper tuning of excitatory input, which amplifies the 
modulation. Therefore, these models cannot be applied 
to the case of a salt-and-pepper structure, as found 
in rodents, with no apparent spatial or feature specific 
connectivity. 

Even in species with orientation maps, there seem to be 
some issues with these models. First, they rely on - and 
predict - a sharpening mechanism of selectivity due to 
tuned recurrent excitation. However, the late (presumably 
recurrent) sharpening of selectivity, which these theories 
predict, has not been observed in experiments (Gillespie 
et al. 2001; Sharon and Grinvald 2002). Also, the ori- 
entation selectivity of neurons seem to be the same as 
their feedforward input, since the preferred orientation 
of neurons does not change with recurrent interactions 
(Gillespie et al. 2001). Rather than sharpening of tuning 
curves, a more plausible function of the recurrent net- 
work is increasing the modulation ratio, by suppressing 
the baseline (Sharon and Grinvald 2002). 

This was indeed the main mechanism of orientation 
selectivity we described in our networks here. As it is 
based on essentially linear processing, our model predicts 
no sharpening of the tuning as a result of recurrent inter- 
actions, but only an increase in modulation depth, not 
affecting the tuning width (Sharon and Grinvald 2002) d . 
Sharpening of tuning curves would only be a consequence 
of the feedforward nonlinearity, reflected in a nonlinear 
transfer function of neurons. As our results do not depend 



on the power-law transfer function of single neurons, our 
model would also work if the operating regime of the cor- 
tex suggested a smaller exponent of the power-law (Xing 
et al. 2011) e . Moreover, as we demonstrated above, this 
mechanism does not have to be accompanied, on average, 
by a large shift between the input and output preferred 
orientations. 

Another consequence of the sharpening theories is the 
emergence of strong pairwise correlations in the network 
(Series et al. 2004). This seems not to be consistent with 
the very low correlations reported in the neocortex (Ecker 
et al. 2010). More specifically, it has recently been shown 
that highly selective neurons in the input layer of monkey 
VI exhibit very low noise correlations 2012. This imposes 
an important constraint on recurrent models which need 
sharper tuning of excitatory input to the neurons as com- 
pared to inhibition, as this sharper tuning leads to higher 
noise correlations in the local network (see Figure five in 
Hansen et al. (2012)). Hence, it might be difficult for these 
models to simultaneously account for both sharp orienta- 
tion selectivity and low pairwise correlations in the input 
layers. 

The mechanism we discussed here, in contrast, does 
not need - and not predict - strong pairwise correla- 
tions in the network. In fact, our mean field analysis was 
even based on the assumption of no correlations in the 
network. As the network operates in the Al state, the 
amount of linear read-out of information from our net- 
works would therefore be several times higher than in 
sharpening theories, comparable to feedforward models 
(Series et al. 2004). 

In comparison to feedforward models, however, our 
analysis suggests that contrast invariant tuning of both 
membrane potentials and spiking activity (Anderson 
et al. 2000) can robustly and reliably emerge through the 
action of a recurrent network. Contrast invariance is a 
critical property of feature selectivity, which ensures reli- 
able and consistent feature detection for a wide range of 
different stimuli. Without a network mechanism of the 
sort described here, neurons would need a specific fine- 
tuning for each contrast, in order to be selective for the 
same feature. The network mechanism proposed here 
provides a generic mechanism to dynamically achieve 
contrast invariance, without the need for feature specific 
wiring, special correlation structure, power-law transfer 
function, contrast-dependent variability, shunting inhibi- 
tion, synaptic depression, adaptation or learning. How- 
ever, it remains to be experimentally verified whether 
intra-cortical recurrent connectivity is indeed necessary 
for contrast invariance, or whether feedforward mecha- 
nisms are enough to account for this phenomenon (Finn 
et al. 2007; Priebe and Ferster 2012; Sadagopan and Ferster 
2012). A crucial experiment would be to test whether the 
tuning of the membrane potential is still contrast invariant 
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if lateral interactions in the cortex are deactivated (Chung 
and Ferster 1998; Ferster et al. 1996; Kara et al. 2002). 

Future studies 

There are, however, several issues which need to be fur- 
ther examined in future works. 

Extending the scope of linear analysis 

First, our analysis is mainly provided to compute the mean 
values of baseline and modulation gains in the network. 
It is therefore necessary to extend the analysis such that 
it accounts for the distribution of these gains. Also, the 
model assumes cosine tuning of inputs (linear tuning), and 
linear network operation (e.g. no rectification in the tun- 
ing curves). It would therefore be revealing to see if, and 
to which degree, the results of the linear analysis hold 
in the presence of nonlinearities reported to exist in the 
biological cortex (Anderson et al. 2000). 

Neuron model 

We used current based LIF neurons with unconstrained 
membrane potentials in our simulations, since this gave 
us the opportunity to perform a theoretical analysis of 
network dynamics. It is however necessary to test to 
which degree the results of our study change by recruit- 
ing a different neuron model. For instance, using an 
alternative neuron model like concuctance-based LIF 
neurons might not allow the distance to threshold of 
the membrane potential to increase unboundedly. This 
was not the case here, as we did not impose any 
minimum bounday condition on our current-based LIF 
neurons. 

Such a difference may change the effective gain of neu- 
rons and, as a result, a different eigenvalue spectrum 
might be obtained. This, in turn, may change network 
dynamics and lead to a qualitatively different behavior 
of the network. It might also have consequences for the 
structure of correlations in the network, and may affect 
the Al state. Our preliminary results indicate that impos- 
ing a lower boundary condition on current based LIF 
neurons can amplify correlations and synchrony in the 
network, to the extent that the network does not operate 
anymore in the Al state. Feature selectivity is still obtained 
and even enhanced in the network. Tuning curves show 
a higher average OSI and maximum firing rates, more 
rectification and reduced tuning widths follow from this, 
while contrast invariance is preserved (not shown). Such 
a scenario should be analyzed in more detail in future 
studies. 

Network connectivity 

It is also important to study the effect of different 
connectivity patterns in the network. Here, we have 



modeled the dominance of inhibition by increasing the 
relative strength of IPSPs. An alternative implementation 
is to increase the density of inhibitory axonal projections, 
which increases the density of connectivity. Dense pattern 
of connectivity has been reported for inhibition in differ- 
ent cortices (Fino and Yuste 2011; Hofer et al. 2011; Packer 
and Yuste 2011), and seems to be a common motif. Such a 
change in network connectivity might have consequences 
for sensory processing. For instance, a decrease in tem- 
poral fluctuations of the local inhibitory population and, 
likewise, in the quenched noise of preferred orientations 
is implied, which can, in turn, affect the tuning of inputs 
and amplification or attenuation of orientation selectivity. 

Also, the model and the analysis provided here should 
be extended to account for networks with spatial struc- 
ture. It should be analytically studied how distance depen- 
dent connectivity, and in particular different connectivity 
profiles for excitation and inhibition, affect the results 
obtained here. It has been shown, for instance, that bal- 
anced networks can show topologically invariant statistics 
(Yger et al. 2011). It would therefore be interesting to 
see if the same analysis also applies to functional prop- 
erties of these networks. Of special interest would be to 
study how the spectrum of the network changes (Voges 
et al. 2011), and to which degree this predicts the opera- 
tion of the network, in particular, how the spatial extent 
of excitation and inhibition affects this behavior. Exper- 
imentally, it has been reported that inhibition is more 
local than excitation in terms of anatomical projections. 
Many theoretical models, however, assume broader inhi- 
bition for convenience. Studying the functional properties 
of networks with realistic patterns of connectivity might 
therefore shed light on this aspect. 

Orientation selectivity and orientation maps 

As we were primarily interested in the emergence of ori- 
entation selectivity in species without orientation maps, 
we studied here random networks with salt-and-pepper 
structure. However, the model could be extended to net- 
works with spatial or functional maps, which imply feature 
specific connections. As opposed to the Mexican hat 
profile assumed in the ring model, if one now assumes 
a more realistic pattern of local inhibition and longer 
range excitation, different dynamic properties might fol- 
low (Hansen et al. 2012; McLaughlin et al. 2000; Pernice 
et al. 2011). The analysis of the new regime of orien- 
tation selectivity therefore calls for a further study. The 
results of this modeling would in turn help identifying 
the basic mechanisms that are responsible for the emer- 
gence of orientation selectivity in different species with 
different structures, and to eventually provide an answer 
to the question whether common design principles exist, 
or whether different strategies have been recruited by 
different species. 
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Methods 

Simulation and analysis of network dynamics 
Neuron model and surrogate spike trains 

We studied networks of leaky integrate-and-fire (LIF) 
neurons. For this spiking neuron model, the sub-threshold 
dynamics of the membrane potential Vt(t) of neuron i is 
described by the leaky-integrator equation 



currents are pulse-like, the dynamic equation for the net- 
work is obtained from 



rV,(t) + Vi(t) = RIi(t). 



(13) 



The current lift) represents the total input to the neuron, 
the integration of which is governed by the leak resistance 
R, and the membrane time constant r = 20 ms. When 
the voltage reaches the threshold at V t h = 20 mV, a spike 
is generated and transmitted to all postsynaptic neurons, 
and the membrane potential is reset to the resting poten- 
tial at Vrj = 0 mV. It remains at this level for short absolute 
refractory period, t le f = 2 ms, during which all synaptic 
currents are shunted. 

To simulate spiking inputs to neurons from outside the 
network (e.g. from the lateral geniculate nucleus, LGN), 
we resorted to the conceptually simpler model of a Pois- 
son process. The associated surrogate spike trains have 
the property that spikes are generated randomly and inde- 
pendently with a prescribed firing rate at each point in 
time. The linear superposition of an arbitrary number of 
Poisson processes (as in the case of multiple afferents) 
is again a Poisson process. The rate of the superposition 
process is exactly the linear sum of the rates of its com- 
ponents, and it can be effectively simulated as a single 
process with high rate. 

Network connectivity and activity dynamics 

The networks considered in this study comprised N = 
12 500 neurons, / = 80% of which were excitatory and 
1 — / = 20% inhibitory. Synaptic connections were 
drawn randomly and independently, such that each neu- 
ron received exactly 1 000 inputs from the excitatory and 
250 from the inhibitory neuron population, respectively. 
This amounted to an overall connectivity of e = 10%, as 
suggested by statistical neuroanatomy of local cortical net- 
works (Braitenberg and Schiiz 1998). The wiring imposed 
in our model was in accordance with Dale's principle, 
i.e. each neuron formed the same type of synapse with all 
its postsynaptic partners, either excitatory or inhibitory 
(Kriener et al. 2008). Self-connections were excluded. 

If t\ is the time of a spike elicited by neuron i, we use 
a Dirac delta-function S(t — tf) to represent it as a time 
dependent signal. The sum Yift) = J2k^^ ~ tnen 
stands for a spike train. The input lift) to each neuron is 
the sum of all excitatory and inhibitory postsynaptic cur- 
rents (PSCs) induced by presynaptic spikes that arrive at 
its dendrite, and the hyperpolarizing currents responsi- 
ble for the reset after each output spike. Assuming that all 



Rlift) = r 



-V Kset W) + WijYjft - D) +J s Xift) 



(14) 



After each spike, the membrane potential was reset to the 
resting potential at OmV, therefore the size of the volt- 
age jump is Vreset = Vth — Vn = Vth- The cross-neuron 
coupling Wij encodes the amplitude of the postsynaptic 
potential (PSP), corresponding to a synapse from neuron /' 
(source) to neuron i (target). A uniform transmission delay 
of D = 1.5 ms was assumed for all recurrent synapses in 
the network. The spike train Xift) stands for the accumu- 
lated external input to neuron i, and the corresponding 
synapses have connection strength of amplitude J s . 

In all the simulations described in our paper, in fact, 
we used stereotyped synaptic transients of finite width, 
instead of normalized impulses, as postsynaptic currents. 
All synaptic kernels had the shape of an alpha-function 
Ia^r-te~^ Ts y n , with a fixed time constant r svn = 0.5 ms, 

~syn ' 

replacing the delta-functions in the spike trains. The peak 
amplitude of the kernel is J a , to which we refer as EPSP 
to denote the strength of post-synaptic potentials. The 
parameter Wij corresponds to the integral of the PSP, 
which is / = er S y n EPSP. 

In this model, keeping all time constants at fixed val- 
ues, the efficacy of a synaptic connection is determined 
by the peak amplitude of the PSP. For any specific net- 
work, we assumed that all recurrent excitatory synapses 
induce excitatory postsynaptic potentials of the same peak 
amplitude, EPSP. The peak amplitudes of inhibitory post- 
synaptic potentials were taken to be a fixed multiple, g, of 
the excitatory ones, such that IPSP = — gEPSP. For all our 
results presented in the main text, individual inhibitory 
couplings were assumed to be much more effective than 
excitatory ones: The excitation-inhibition ratio was fixed 
at g = 8. As a consequence, recurrent connectivity in 
our networks was characterized by a net surplus of inhi- 
bition, since the small number of inhibitory neurons was 
over-compensated by the strength of individual inhibitory 
couplings. Fixing the balance between recurrent excita- 
tion and inhibition in this way is an important concept in 
models of cortical dynamics, although measurements in 
real brains are difficult (see e.g. Okun and Lampl (2008)). 

In different simulations, we used excitatory synapses 
with an EPSP amplitude in the range between 0 mV and 
1.0 mV. Accordingly, inhibitory synapses had IPSP ampli- 
tudes between OmV and — 8.0 mV. All external inputs 
in our simulations were excitatory, and the amplitude of 
their synapses, EPSPff w , was fixed at 0.1 mV throughout all 
simulations. 
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This configuration of parameters, combined with a sta- 
tionary driving input to each neuron in the network, 
was previously shown to induce relatively low rates in 
all neurons, while spike trains are irregular, and pair- 
wise correlations remain weak (Brunei 2000). These prop- 
erties are known to be a result of complex recurrent 
network dynamics, and not a consequence of random 
inputs (e.g. Poisson spike trains) that drive the network 
(Kriener et al. 2008; van Vreeswijk and Sompolinsky 
1996). Inhibitory feedback actively decorrelates the net- 
work activity (Pernice et al. 2011; Renart et al. 2010; 
Tetzlaff et al. 2012). The resulting states of network 
dynamics are dubbed asynchronous-irregular, Al, and 
they are thought to closely resemble the dynamic states 
of neocortical networks recorded with extracellular elec- 
trodes (Ecker et al. 2010). 

In this parameter setting, the degree of recurrence in 
any specific network is essentially determined by the 
amplitude of excitatory postsynaptic potentials, EPSP, 
of the recurrent connections. Recurrence can be effec- 
tively quantified by the spectral radius of the connectivity 
matrix W, which scales linearly with the EPSP amplitude. 
This fact is explained in more detail below. 

Neuronal tuning 

To explore the tuning curves of neurons in a network, we 
simulated their responses to stimuli with different orien- 
tations. Beyond excitatory and inhibitory input from the 
recurrent network, each neuron received extra excitatory 
input, the firing rate of which exhibited a slight depen- 
dence on stimulus orientation. This external input can 
be conceived as the overall effect of stimulation, and it 
includes inputs from LGN, and possibly afferents from 
other, non-local cortical sources. In our simulations, the 
input was implemented as a homogeneous Poisson pro- 
cess, with an average firing rate, s, depending on the 
stimulus orientation 6 according to 

s(6) = s B [l + mcos(2(e -6*))]. (15) 

The baseline sg is the mean level of input across all ori- 
entations. We used a logarithmic relation between input 
contrast C and input baseline, sb oc log 10 (l + 100C), as 
a practical way to specify the input intensity, inspired by 
biology. 

In all our simulations, the relative amplitude, m of the 
stimulus dependent modulation was fixed to a fraction 
of 10% of the baseline level, corresponding to setting 
m = 0.1. The parameter 0* signified the stimulus ori- 
entation at which the neuron received its maximal input, 
•Smax = (1 + m)SB- It represented the initial preferred ori- 
entation, Input PO, of the neuron, a parameter that was 
randomly and independently assigned to each neuron in 
the population. 



To measure the output tuning curves in numerical sim- 
ulations, we stimulated the networks for 12 different stim- 
ulus orientations, covering the full range between 0° and 
180° in steps of 15°. Each simulation was run for 6.3 s, 
using a simulation time step of 0.1 ms. In order to include 
only steady state activity into our analysis, and to avoid 
onset transients, the first 300 ms in each simulation were 
discarded. The output tuning curve of any neuron in the 
network was obtained in terms of its average firing rate r 
in response to each stimulus orientation 9, and normally 
plotted as a curve r(Q). An output tuning curve would 
be termed contrast invariant, if its overall shape does not 
depend on the contrast, C, of the stimulus. 

To explore the interaction between feedforward and 
recurrent connectivity on orientation selectivity, we sys- 
tematically changed two parameters in our networks: The 
mean input firing rate, sb, and the EPSP amplitude as 
a measure for the strength of recurrent coupling. We 
changed these two parameters in a network, while fix- 
ing all other parameters, including the network topology 
given by a specific realization of the random synaptic con- 
nectivity, W, the inhibition-excitation ratio, g, and the 
input modulation ratio, m. We used three different val- 
ues for the baseline intensity sb = 12000, 16000, and 
20000spikes/s. This is corresponding to low, medium, 
and high contrast, C « 9%, 39%, and 99%, respectively. 

Free membrane potential 

Contrast invariance of the membrane potential tuning is, 
in the case of a tuned spike response, compromised by the 
reset mechanism in our neuron model: After each spike, 
the membrane potential is reset to its resting value, which 
exerts a negative contribution to the membrane potential, 
which effectively imposes the opposite tuning as com- 
pared to the output spiking. As a result, when the neuron 
fires more (higher contrasts), it inevitably attains a more 
negative membrane potential. To correct for this phe- 
nomenon, we add this negative contribution back to the 
membrane potential (Vf m = V+rrV Kset ) or, equivalently, 
keep the neuron from spiking by raising its threshold 
to very high levels. In membrane potential recordings 
from real neurons, tt is also common to correct for this 
phenomenon by cutting out the spikes including their 
aftereffects. If used with care, this procedure is essentially 
equivalent to the correction we applied here. 

Numerical methods and simulation software 

The implementation of the LIF model employed in the 
present study is based on a numerical method known 
as "exact integration" (Diesmann et al. 2001; Rotter and 
Diesmann 1999). Numerical simulations of all networks 
were performed using the neuronal simulation environ- 
ment NEST (Gewaltig and Diesmann 2007). This tool 
has been developed to support the reliable, precise and 
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performant numerical simulations of networks of spiking 



neurons. 



Data analysis 
Orientation selectivity 

To quantify orientation selectivity, we computed the Pre- 
ferred Orientation, PO, and the Orientation Selectivity 
Index, OSI, of each neuron from its respective tuning 
curve, r(6), obtained in numerical simulations. To this 
end, we first computed the circular mean (Batschelet 
1981) of the firing rate measured at each orientation, 
which we call the Orientation Selectivity Vector, OSV, 



OSV: 



J^o r{9) exp(27r/67180°) 



(16) 



The PO, 9*, was then extracted as the angle, arg(OSV), 
of the OSV. Its length, |OSV|, yielded a measure for the 
degree of orientation selectivity, OSI (Ringach et al. 2002). 
For a highly selective neuron, which is only active for one 
orientation, and remains silent for all other orientations, 
the OSI would be 1; for a completely unselective neuron 
responding with the same firing rate for all orientations, 
this measure returns 0. 

For better comparison with the experimental literature 
(see e.g. Niell and Stryker (2008)), an alternative measure 
of orientation selectivity has also been computed for the 
tuning curves obtained in our simulations. It is given by 



OSI* = 



f pref I" orth 
? pref H~ ? orth 



(17) 



where r pre f = r(9*) is the firing rate at the preferred stim- 
ulus orientation, 9*, and the r or th = r(G* +90°) is the firing 
rate for the orientation orthogonal to it. 

Since the output preferred (and hence the orthogonal) 
orientations of a neuron are computed from Eq. (16), we 
need to interpolate between the data points of a tuning 
curve to obtain r pre f and r ort h. To do this, we fit a cosine 
function to the tuning curve sampled at 12 equidistant 
orientations, employing a nonlinear least squares method. 
Then, the cosine fit of the tuning curve is evaluated at PO 
and PO + 90° to obtain r pre f and r or th> respectively. Neg- 
ative numbers, whenever occurring, were replaced by 0. 
This is similar to what experimentalists typically do (see 
e.g. Niell and Stryker (2008), although by fitting different 
functions, like a Gaussian or a von Mises probability den- 
sity; also see Hansel and van Vreeswijk (2012)), therefore 
allows us to compare our distributions to their reported 
results f . 

For a perfect cosine tuning curve according to Eq. (15), 
we obtain 



OSI 



m 



and OSI* = m 



(18) 



irrespective of the baseline firing rate sg, and of the pre- 
ferred orientation 9*. Thus, in the case of our input tuning 
curves with m = 0.1, we have OSI = 0.05 and OSI* = 0.1, 
respectively. 

Baseline and modulation gain 

To quantify the processing of baseline and modulation in 
a specific network, we compute the mean baseline and 
mean modulation gains over all neurons. To obtain this, 
we first compute baseline and modulation gain for indi- 
vidual tuning curves, r„(9), as follows. The baseline is 
obtained by averaging the tuning curve over all orienta- 
tions 



1 K 
k=l 



(19) 



where r n {9k) is the firing rate of «-th neuron in the net- 
work in response to a stimulus with the A~-th orientation, 
and K = 12 is the number of different orientations 
considered in the simulation. We refer to this as the F0 
component of tuning curves. The modulation is conve- 
niently obtained from the absolute value of the second 
Fourier component of the tuning curve 



^ = I 5>«(%)exp(-27U0/18O°)|. 



(20) 



k=l 



We refer to this as the F2 component of a tuning curve. 

The baseline gain and the modulation gain are then 
defined as the output value divided by the input value, 
respectively 



Yr. 



= t*/4 and y? = «. 



(21) 



The input baseline and modulation, and s^, are 
obtained from input tuning curves in the same fashion as 

and were obtained from the output tuning curves, 
respectively (Eq. (19) and (20), respectively). The corre- 
sponding mean values, yb and ym> for any network are 
the average over all individual gains. Finally, each gain is 
normalized by its value for a network with no recurrence, 

and Ym' which are obtained from the simulation of a 
network with no recurrence, EPSP = 0 mV 



, , norm Y B 

Yb =~o 
Yb 



and 



norm YM 



(22) 



PO scatter 

The transformation of preferred orientation induced by 
a recurrent network is visualized by a scatter plot show- 
ing Output PO vs. Input PO for all neurons (Figure 12). 
Weakly recurrent networks essentially preserve the pre- 
ferred orientations of the input to each neuron, leading 
to scatter plots centered about the diagonal. For networks 
with increased recurrence, the output PO deviates from 
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the input PO, and off-diagonal elements occur more fre- 
quently. To quantify the deviation, we first compute the 
difference, APO„ = OutputPO H - InputPO„, for each 
neuron. Observe that orientation should be taken mod- 
ulo 180° and APO„ = 90° represents the largest possible 
difference between input and output PO. 

As a numerical measure for the total degree of PO scat- 
ter in a network, we computed the Scatter Degree Index, 
SDI (Figure 12B). Its definition is based on the circular 
mean 



1 N 

R A = — ^ exp(2ni APO„/180°), 



(23) 



n=\ 



where N is the number of neurons in the network. The 
SDI is then given by the angular deviation (Batschelet 
1981), which can be computed from the length \R&\ 
according to 



90° , 

SDI = — VAl 

71 



\Rk\). 



(24) 



Note that as APO spans the half-cricle, i.e. the range 
[0°, 180°], we have taken half the resultant angle as the 
SDI. If all Output POs are exactly the same as Input POs, 
SDI returns zero; the maximum scatter from the Input 
POs corresponds to a uniform distribution of APO, for 
which SDI returns « 40.5°. 

Theoretical analysis 

Firing rates and membrane potential statistics 
Mean firing rates in recurrent networks 

The tuning curves considered in this work reflect time- 
averaged firing rates of neurons in a recurrent network. 
From our numerical simulations, it became clear that the 
time averaged membrane potential is indicative of the 
operating point of the network with regard to the tuning 
properties of its neurons (see Section Results). There- 
fore, we begin our analysis by considering time averaged 
equations. 

Assuming stationarity, we form temporal averages ( . ) 
of all dynamic variables that occur in Eq. (13) and (14). 
Since there can be no drift of the time averaged membrane 
potential in this case, we have (Vf) = 0. We write v,- = 
{Vi) for the time averaged membrane potential, r; = (Yi) 
for the mean firing rate of neuron i in the network, and 
Si = (Xj) for the mean firing rate of its external input, 
respectively. We obtain an equation that relates the sta- 
tionary firing rates of all neurons in the network with their 
mean membrane potentials (Hansel and van Vreeswijk 
2002; Kriener et al. 2008) 



Vtm + ^Wijrj+JsSi 



Observe that transmission delays do not matter for tem- 
poral averages, and that the above equation holds for 
networks of LIF neurons with arbitrary connectivity. 

From now on, we rescale Eq. (25) such that all fir- 
ing rates are expressed in units of 1/r, and all voltages 
are given in units of Vth- For a network of N neurons, 
the recurrent synaptic connectivity is encoded by a fixed 
N x N coupling matrix W = (Wij). The external inputs, 
the firing rates and the membrane potentials of all neurons 
are represented by the AT-dimensional vectors s = (s;), 
r = (ri) and v = (v;), respectively. The time averaged 
equation above then reads, in matrix-vector notation 



-f+ Wr+] s t 



(26) 



(25) 



Solving for the vector r of recurrent firing rates, we obtain 

r = A(J s s - v), where A = (I - W)~ l . (27) 

We assume that the matrix 1 — W is always invertible, 
with inverse A. If two out of the three variables s, r and 
v are known, the third one can then be computed in a 
straightforward fashion. 

Eigenvalue spectrum of homogeneous random networks 

We now specifically consider a recurrent network of exci- 
tatory and inhibitory neurons, as discussed above. It is 
assumed that Ne = JN neurons are excitatory, forming 
synapses of uniform strength / with their postsynaptic tar- 
gets. The remaining neurons Nj = N — Ne = (1 —f)N are 
inhibitory, forming synapses of uniform strength —gj. The 
factor g > 0 describes the relative strength of inhibitory 
synapses. We refer to the network as being "inhibition 
dominated", if the lower number of inhibitory neurons 
is compensated for by stronger inhibitory weights. In the 
case considered here, this amounts to the condition g > 
Ne/Nj =//(l -/) (Brunei 2000). 

The connectivity of the network is set to e, such that 
each neuron receives input from exactly Ne€ excitatory 
neurons and Nj€ inhibitory neurons. The presynaptic 
sources are randomly selected from the available pool, 
multiple synaptic contacts are excluded. The graph under- 
lying such a network is a specific type of random graph 
(Bollobas 2001; Erdos and Renyi 1959). The connectivity 
matrix W is a random matrix with two types of entries, 
organized in homogeneous columns. The entries in pos- 
itive columns of this matrix, corresponding to excitatory 
neurons, have a mean of i]e = Je and a variance of 
erj = / 2 e(l — e), whereas the entries in negative columns, 
corresponding to inhibitory neurons, have a mean of 
rji = —gfe and a variance of erf = g 2 J 2 €(l — e), respec- 
tively. The matrix W has an eigenvalue spectrum with 
two components that are, for large and not too sparse 
networks, described as follows (Rajan and Abbott 2006): 
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There is one exceptional eigenvalue, proportional to the 
mean recurrent input to each neuron 



X 0 = Neve + N m = JNe [f - g(l -/)] . 



(28) 



It belongs to uniform eigenvectors with all components 
being equal. They represent a 1 -dimensional subspace, 
spanned by the uniform vector u = In a 

balanced random network, we have An < 0. The bulk 
spectrum A = [k\, . . . , A.jv-i} covers a circular region 
in the complex plane, centered at the origin. Its radius p 
satisfies 

p 2 = N E a 2 + Nrf = J 2 N€(1 - e) [f+g 2 (l -/)] . 

(29) 

The density of eigenvalues within the circle is in gen- 
eral non-uniform, and it can be approximated by a density 
derived in Rajan and Abbott (2006). 

Eq. (26), which relates input, output and membrane 
potentials under stationary conditions, has an effective 
coefficient matrix W— 1. Its eigenvalue spectrum consists 
of numbers X — 1, where X is from the spectrum of W. 
Likewise, the eigenvalues of (1 — W)~ x are (1 — X)~ l . This 
can either be derived directly, or it can be implied by the 
spectral mapping theorem (Higham 2008). The associated 
eigenvectors are the same in each case. 

Self-consistent firing rates in homogeneous random networks 

Under the same conditions on homogeneity as made 
above, explicit solutions for the response rates, and for the 
mean membrane potentials, can be obtained by resort- 
ing to additional constraining assumptions. Specifically, 
one can analytically describe the response statistics of a 
leaky integrate-and-fire neuron, which is driven by ran- 
domly fluctuating input. If inputs are uncorrelated, and 
synaptic couplings are weak, the lumped synaptic input 
current may be approximated by a Gaussian white noise 
with appropriate parameters \i and a 2 



RI(t) = n + Oyfrriit) 



(30) 



where rj{i) is a stationary Gaussian white noise with zero 
mean and unit power spectral density. Assuming sta- 
tionarity and a fixed voltage threshold, the associated 
first-passage time problem can in fact be solved: The 
membrane potential dynamics of the neuron can be con- 
ceived as a diffusion process, and the time evolution of 
the membrane potential distribution is given by a Fokker- 
Planck equation with specific boundary conditions. Its 
solution yields explicit expressions for the moments of the 
inter-spike interval distribution (Ricciardi 1977; Siegert 
1951). In particular, the mean response rate of the neuron, 



r, in terms of its input statistics 
r =f(fi,a) 

t re f + I e" (1 + erf(w)) du 



(31) 



with V t h = (Vth - M)/ cr and V 0 = (Vrj - ^)/cr. 

Employing a mean field ansatz, the above theory can be 
applied to networks of identical pulse-coupled LIF neu- 
rons, randomly connected with homogeneous in-degrees, 
and driven by external excitatory input of the same 
strength. Under these circumstances, all neurons exhibit 
the same mean firing rate, which can be determined by 
a straight-forward self-consistency argument (Amit and 
Brunei 1997; Brunei 2000): The firing rate r is a func- 
tion of the first two moments of the input fluctuations, /x 
and a 2 , as described by Eq. (31). Both parameters are, in 
turn, functions of the firing rate r. This leads to a fixed 
point equation, the root of which can be found numeri- 
cally. Here we employed Newton's method, verifying the 
convergence of the iteration by appropriate means. 

For networks of the type described here, we have specif- 
ically 



2 



T[J s s+JrNe(f-g(l -/))], 
T[J 2 s + J 2 rNe(f + /(!-/))], 



(32) 



where s is the input (stimulus) firing rate, and r is the mean 
response rate of all neurons in the network, respectively. 
Here, J s denotes the EPSP amplitude of external inputs, 
and / denotes the amplitude of recurrent EPSPs. The 
inhibition-excitation ratio g has been introduced above. 
The remaining structural parameters are the number of 
neurons in the network, N, the connection probability, 
e, and the fraction / of neurons in the network that are 
excitatory, implying that a fraction 1 — / is inhibitory. 

Correction for a-synapses 

The treatment described above is only approximating the 
networks considered in numerical simulations, since we 
chose biologically more realistic LIF neurons with alpha- 
synapses. In order to make use of the same analytical 
framework as just described, we made the simplifying 
assumption that all the presynaptic current is delivered 
immediately, and that the input current to each neuron is 
still white. We therefore need to obtain the effective values 
for mean and variance. 

To obtain the effective value of the mean, we match the 
area under the PSC kernel of a-shape with a correspond- 
ing S-synapse 



f 

Jo 



t 

■"syn 



" f /%" dt = l. 



(33) 
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The actual a-synapse with a peak amplitude } a would then 
be matched to a <5-PSC as follows 



f 

Jo 



J a te t/Ts y n dt = 



(/a^^syn) 



L syn 



/ —e~ t/T ^ n dt. 

J 0 T syn 



(34) 



Therefore, we choose the value / = / a er S y n as the effec- 
tive value for the mean input. This is equivalent to the 
integral under the PSC, i.e. the total amount of current 
that is delivered by an alpha synapse with peak J a . 

The effective value of the variance can be obtained in 
the same fashion by matching the integral of the squared 
PSC, of the a-PSC with the S-PSC 



L 



'•syn 



[te 



-t/x, 



vm f dt : 



1, 



(35) 



and 



f \J a —te 

JO L r syn 
JO 7 svn 



-t/x S) 



n2 



-t/x s) 



dt 



| dt. 



(36) 



This suggests J z = \jle 2 r syn . 



Transformation of tuning by a recurrent network 
Baseline and modulation 

We now consider the case of tuned input to a recurrent 
network. The input S; to each neuron in the network, as 
well as its firing rate response r; and its membrane poten- 
tial Vi, may depend on a given feature 4> e F of a sensory 
stimulus. The functions s;(</>)> r i(<P) an d V/(0) will be called 
the input, output and membrane potential tuning curves 
of neuron i, respectively. Let now the input to each neuron 
in a recurrent network be tuned, i.e. s = s(<fi). After relax- 
ation to equilibrium, the output of the recurrent network 
is given by Eq. (27), and the tuning curves are obtained by 



H4>) = (I - wr 1 [/ s s(0) - v(4>)] 
= A[j s s(4>)-v(4>)]. 



(37) 



Assume now that a stimulus ensemble has been fixed 
for an experiment, think of a uniform distribution for the 
orientation of a stimulus offered, for example. Mathemati- 
cally, this is described by a suitable probability distribution 
on the set of possible values for the feature (/>, the stimu- 
lus ensemble. Thereby, any stimulus dependent quantity x 
(like s, r and v) turns into a real- valued random variable. Its 
expected value E[x] then corresponds to the component 
of the tuning curve x(4>) that is common to all parame- 
ter values, the baseline of the tuning curve. Note that this 



concept depends on the stimulus ensemble, and it sug- 
gests the following further terminology: 



baseline : 
modulation : 



xb - 

Xm 



E[x] 



■E[x] 



(38) 
(39) 



Evidently, the decomposition x = xb+xm is fully specified 
by these settings. Moreover, by linearity of E[.] Eq. (27) 
implies 



h =A(J s s B - v B ), 
and, therefore, 

m = r - r B = A [j s (s - s B ) - (v - v B )] 
= A(J s s M - v M )- 



(40) 



(41) 



This means that the recurrent network defined by the 
coupling matrix W, and the matrix^, derived from it, pro- 
cesses baseline and modulation components separately 
and independently, with no cross-talk involved. In other 
words, pure modulation input will not attain any baseline 
through network processing, and vice versa. This is exactly 
the meaning of Figure 4B. 

Note, however, that, as the mean membrane potential v 
actually depends on the input s in a highly nonlinear fash- 
ion, the above equations determine the network response 
only implicitly. Moreover, for the processing to be inde- 
pendent, it is necessary that vb and v M depends only on sb 
and smi respectively, with no cross-talk. For the baseline 
firing rates, this implies that v B is not affected by the mod- 
ulation in the input. As baseline firing rates are the same 
in our homogeneous networks, the mean vb (over neurons 
in the network) should therefore be more or less constant 
in one experiment with fixed sb- We have checked this 
numerically in our simulations by plotting the standard 
deviation (over neurons) of the mean membrane potential 
(over time and over orientation) for 24 sampled excitatory 
and inhibitory neurons (Figure 16). The variance is indeed 
much smaller than vm, the modulation of the membrane 
potential due to modulation in the input, sm, and this is 
consistent for all recurrent regimes. 

Baseline attenuation by inhibition dominated random 
networks 

We will now consider the case of homogeneous properties 
of all inputs to the network. Specifically, we assume that 
the input baseline is a uniform vector 



E[3] = sb ~ u. 



(42) 



We now consider an inhibition dominated random net- 
work as described in Section Eigenvalue spectrum of 
homogeneous random networks, where all neurons have 
identical parameters, and the connectivity matrix is statis- 
tically homogeneous. As we have already discussed, under 
the assumption of orthogonality, the firing rate responses 
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of all neurons to a homogeneous stimulus, and the cor- 
responding mean membrane potentials, are all identical 

r B ~ u and vg ~ u. (43) 

In other words, homogeneous vectors, like baseline input, 
are eigenvectors of the coupling matrix W. They are 
also eigenvectors of the matrix A = (1 — W)~ l which 
determines the stationary firing rates. The corresponding 
eigenvalue An of the matrix W is in fact negative, and the 
corresponding eigenvalue of A is 1/(1 — An) is positive, 
but much smaller than 1, since |An| is typically large (of 
order N/e). The overall amplification or attenuation of the 
baseline is given by 

h = A(Jsh - vb) = - 1 - (Jsh - v B ). (44) 

1 — AO 

In the case of uniform input sg an explicit solution of 
the mean firing rate rg can be obtained by the mean field 
approximation, as described in Section Self-consistent 
firing rates in homogeneous random networks. The mean 
membrane potential vg is then determined by Eq. (26). We 
will now discuss how the modulation part of tuned inputs 
can be approximated. 

Modulation of the firing rate and of the membrane potential 

We used the theory described in the previous sections to 
approximately determine the response of our networks, 
when they are processing non-uniform inputs, tuned to 
some stimulus feature. Specifically, the mean (over the 
network) modulation component (F2 component) of the 
output tuning curves in response to a tuned input can 
be obtained approximately. Under the assumption of sta- 
bility (see Section Operating regime of network), and for 
the tuned inputs considered here, it seems justified to 
start with the approximation of "perfect balance" of the 
recurrent modulations: 

Wr M « 0. (45) 

Although this approximation is not strictly true, it holds 
on average: As the result of numerical simulation in 
Figure 9B demonstrated, Wsm had a narrow distribution 
around zero. Similar distribution is expected for W?m> as 
?m can be expanded in terms of powers of sm (Eq. (6)) 
under the assumption of stability. 

In terms of diagrammatic illustration of Figure 4B, this 
is equivalent to Pm ^ 0. This in turn implies that 
the net input from the recurrent network is on average 
untuned. For EPSP = 0.2 mV, this input tuning is shown in 
Figure 10. Although the tuning of recurrent input is (com- 
pared to feedforward tuning) not negligible for all neu- 
rons, it holds on average, such that average tuning curves 
have the same shape as the input tuning (Figure 10B). We 
therefore use this approximation to compute the mean 
modulation gain of the network. 



Under this assumption, the computation of modulation 
rate vector, tmi is simplified to 

m « JsSm - v M - (46) 

The linear mixture of tuning curves described by Eq. (37) 
is reduced to an amplification or attenuation of the respec- 
tive input tuning curves. 

Since vm depends nonlinearly on input parameters, 
we again need compute the self-consistent firing rates 
employing mean field theory. Mean /x and variance a 2 
of the input current analogous to Eq. (32), however, are 
now computed by approximating the recurrent firing rate 
by the baseline firing rate, rg, which is the same for all 
neurons in the network, as discussed above. The external 
input to this specific neuron, in contrast, experiences a 
feature specific modulation 5 = sg + sm- As a result, we let 

li = T[Js(sB+s M )+JrBNe(f-g(l-f)]), 

a 2 = r[jf(sg+SM) +J 2 rgNe(f+g 2 (l -/))]. (47) 

The parameters J s , J, g, N, e and / are the same as above. 

The resultant firing rate of the neuron according to 
Eq. (31) now differs from its baseline firing rate. The dif- 
ference is, in fact, a good estimate for the modulation in 
the output firing rate of this particular neuron. Due to 
our general homogeneity assumptions, all neurons in the 
network will have the same output modulation, notwith- 
standing the fact that they all have different preferred 
orientations. Note that, to be consistent, the correction of 
/ due to the refractory period should be performed based 
on the modulated firing rate. This becomes specifically 
important for low recurrences, where the modulation rate 
is higher and, as a result, the effect of shunting of input 
due to refractory period (Jrt re {) becomes more prominent. 

If the firing rates have been determined self- 
consistently, Eq. (26) yields the corresponding membrane 
potentials directly. This is true for both the baseline and 
for the modulation, which are defined in the obvious way 
also for the membrane potential. 

Operating regime of network 

The modulation gains can also be computed by linearizing 
the dynamics around the baseline operating regime of the 
network. 

Our results revealed that the mean modulation gain, 
Ym> in the network depends on the mean distance of the 
membrane potential from the threshold (Figure 18). Sub- 
threshold modulations (with regard to the mean-driven 
threshold) were capable of eliciting output firing activity, 
and the input-output relationship (the gain) was inversely 
proportional to the distance to threshold. One way to 
interpret these results is to summarize them in terms of 
the mean and standard deviation of the input that a neu- 
ron receives on average from the network, in the baseline 
state. To this, and alternatively to the mean and standard 
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deviation of membrane potential, which is uniquely deter- 
mined by the input, we refer as the operating point of the 
network. 

The effect of mean, /xg, and standard deviation of input, 
ag, on the modulation gain can be described, respec- 
tively, as shifting the mean membrane potential (and 
hence determining the mean distance to threshold), and 
smoothing (linearizing) the /-/ curve (Anderson et al. 
2000; Miller and Troyer 2002). The linearized gains can 
then be obtained by perturbing the input around the base- 
line state, as it was described in Section Tuning and 
invariance of membrane potentials and Linear tuning in 
recurrent networks. This total embedded gain of the neu- 
ron in response to this perturbation determines the effec- 
tive coupling strength and, as Figure 19 demonstrated, 
predicts the mean modulation gain in these networks 
quite well. 

The embedded gain modulates both feedforward and 
recurrent couplings. The fact that recurrent connections 
are now effectively weighted by these linear gains might 
suggest an explanation why the network exhibits stable 
activity even for highly recurrent regimes. As Figure 6A 
showed, if the spectrum of W is computed from the 
weight matrix normalized by V t h, the radius of the bulk 
of eigenvalues, p, would be larger than one already for an 
intermediate recurrent regime (EPSP = 0.2, / = 0.27 mV, 
and p = 1.68 in this example). If one now computes 
the normalized radius by weighting the coupling strength 
according to linear gains (i.e. Wy = f/, instead of //Vth), 
the new normalized radius, p norm , is not unstable any- 
more (p norm « 0.87 in this example). This coincides with 
our observations of the numerical simulations. 

This enhanced stability has indeed been demonstrated 
to be the case for all networks we have studied here 
(Figure 19, inset). If one now add to this that £ is inversely 
proportional to distance to threshold, it follows that the 
network dynamically settles in a regime of operation 
which stabilizes the bulk of eigenvalues. This is due to the 
fact that in inhibition dominated networks, increasing the 
recurrent coupling also increases the negative feedback 
within the network, which results in more hyperpolarized 
average membrane potentials of the neurons. 8 This in turn 
leads to a smaller relative contribution of each spike from 
a presynaptic source to the firing activity of the postsynap- 
tic neuron, since the distance to threshold has effectively 
increased. The overall increase or decrease in the effec- 
tive gain, f /, depends on how exactly the mean membrane 
potential, v, is affected by / and how f is in turn depending 
on v. 

Linear tuning in the strongly recurrent regime 

For networks with weak to intermediate recurrence, the 
assumption of "perfect balance" allowed a rather accu- 
rate prediction of the gain for the tuned part of network 



activation ("modulation"). This assumption, however, fails 
in the case of strongly recurrent networks. Under the con- 
straints of "linear tuning" some aspects of the problem can 
be nevertheless treated. 

Vectorial features and linear tuning 

We now consider stimulus features that can be repre- 
sented by vectors 0 in F C W for some n > 1. 

Example 1. The direction of a moving light dot stimu- 
lus in the visual field is represented by an angle in [ 0, 2jt) 
or, alternatively, by a vector in F = S l C R 2 , the 1- 
dimensional sphere. The speed of the movement can be 
considered simultaneously with its direction, if encoded 
by the length of the vector. Any vector in R 2 is then 
corresponding to a valid stimulus. 

Example 2. The orientation of a moving grating in the 
visual field corresponds to vectors in F = S 1 C R 2 via the 
bijective mapping 

[0, it) -> S 1 , 8 k> (cos(20), sin(20)) . (48) 

The factor 2 in the argument of the cosine and the sine 
function makes sure that a rotation of the grating by tt 
is mapped to the initial orientation again. Note that this 
parametrization already implies that "cross-orientation 
suppression" follows from the obvious relation &>(0 + 
jt/2) = -co(<t>). 

Example 3. A stimulus for studying color vision is rep- 
resented by the activation profiles of the different types 
of receptor cells in the retina, distinguished by their spe- 
cific light absorption spectrum. For example, trichromacy 
in humans and closely related monkeys involves the dif- 
ferential activation of the three different types of cones 
(S,M,L). This leads, in a natural way, to a representation 
of a color stimulus in terms of a vector in R 3 . 

A simple but relevant model of specific tuning curves is 
linear tuning 

Tiih = ft + (49) 

The parameters ty* > 0 and 0* € F are fixed and spe- 
cific for each neuron: 0* is the baseline rate in absence 
of stimulation, the vector 0* is the preferred feature. In 
fact, stimulating with 0 = 0* produces the highest, and 
stimulating with 0 = —0* the lowest firing rates. More 
generally, if § denotes the angle between the vectors 0* 
and 0, linear tuning is equivalent to cosine tuning 

(0*,0> = ||0*||||0||cos(^). (50) 
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The length of the vector that represents the preferred 
feature ||0* || is the tuning strength, it satisfies 

Wl\\=m&x(if,$)/W\\ (51) 

where ||0|| is the stimulus strength. To ensure that the 
firing rate 7/(0) remains positive 

O<min7'/(0) = 0f-||0f||||0|| (52) 

<t> 

the strength of the admitted stimuli must be limited, so 
we admit only stimuli that are weak enough such that lin- 
earity of the tuning and positivity of firing rates remain 
compatible 

F={0eR":||0||<^*/||^||}. (53) 
Linear tuning in recurrent networks 

Assume now that the individual inputs s/ are tuned with 
respect to an w-dimensional feature 0. The same stimulus 
is "seen" by all neurons, but each neuron responds with its 
private tuning curve 

S/(0) = 7/(0). (54) 

As described by Eq. (37), the responses of neurons in a 
recurrent network have tuning curves that are, in general, 
linear sums of the tuning curves of the input channels. If 
the recurrent interactions are strong, many input channels 
contribute indirectly to the tuning of every output chan- 
nel, recruiting multi-synaptic pathways. Assume now that 
all inputs S; are linearly tuned to the stimulus 0 according 
to 

s/(0) = 0f + <0f,0) (55) 

for parameters 0* and 0*. Then, exploiting the two-fold 
linearity, we obtain 

r(0)=l[/ s s(0)-v(0)] 

= A[j s (jr* + $>*$) -v(&] 

= AJ s jr* + i/ s O*0 - Av{h (56) 

where 0* is the vector of baseline activities and <£>* is 
a matrix the rows of which are given by the transposed 
preferred features (4>*) T ■ Therefore, apart from nonlinear 
distortions induced by nonzero mean membrane poten- 
tials, all neurons in the recurrent network are again lin- 
early tuned, with baselines given by the components of the 
vector A/ S \j/*, and preferred features encoded by the rows 
of the matrix AJ S Q>*. 

Mixing of preferred features in random recurrent networks 

Because linear tuning curves are linearly transformed 
according to Eq. (56), we can actually compute the recur- 
rent preferred features that result from this transforma- 
tion. Note, however, that the actual tuning curves will, 



in general, be contaminated by nonlinear distortions by 
Av(<p) that are not reflected by the linear mix AJ S &* of 
preferred features of the inputs to the network. To keep 
the present discussion simple, we ignore this complication 
here. 

A network with zero recurrent interaction would be 
described by the matrix / s l. Therefore, the perturbation 
of each neuron's private preferred input feature resulting 
from recurrent network action is given by AJ s <i>* — / s <t>* = 
RJ S Q>*, where 

R=J s A-J s l=J s W{l-W)-\ (57) 

If we can assume that the preferred features <t>* of inputs 
are all chosen independently from some common distri- 
bution, the sums of preferred vectors that result from the 
action of the matrix R will, according to the Central Limit 
Theorem, be normally distributed vectors in R". In that 
case, it suffices to compute the covariance matrix C p tb 
of the perturbations of the preferred features of all out- 
put tuning curves, and to compare it with the covariance 
matrix Q n of the inputs 



Cin — 


V 


(3>*) r O* 


Cptb 


1 

N — 


-(RJ s $*) T (RJs®*) 




J 2 
N- 


(®*) T (R T R)$>*. 



Specifically, if the distribution of input features is isotropic 
with covariance Q n = cr^l, the scalar af n represents the 
mean squared tuning strength of all inputs. By means of 
the matrix R, the distribution of output features will then 
be perturbed by a component that is normally distributed 
with isotropic covariance C pt b = o 2 ih \ where 

<b = / S Vn With P 2 = J2 R l ( 59 ) 

/' 

For a random network as described above, the factor 
is the same for all rows /', and it describes now the 
mean attenuation/amplification of the tuning strength 
performed by the recurrent network. 

As an example, we compute here the distribution of 
tuning strengths ||0*|| of neurons that would result in a 
strongly recurrent network, where the perturbation dom- 
inates the result. It is given by the probability density 

<l> n {x;o 2 ) = r tn {x)-f n {x;o 2 ) (60) 

where r] n (x) is the surface of the (w — l)-sphere with radius 
x, and 0„(#; a 2 ) is the probability density function of the 
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w-dimensional normal distribution with zero mean and 
isotropic covariance matrix er 2 l. We find 



r "/2 



<p„(x) = 2 



„«-i 



T(«/2) (2n) n l 2 a n 
A x n ~ l e~ x /(2<T \ 



(61) 

T(n/2)a" 

In the case n = 2, which is particularly relevant for the 
application described in this paper, we get 



b 2 (x;a ) = -^e 



-x 1 /(2a 1 ) 



(62) 



which is a Weibull distribution with shape parameter 2 
and scale parameter ~j2o. 

For the scenario, where feedforward and recurrent com- 
ponents of tuning are mixed, one obtains an interpolation 
between "purely feedforward" and "purely feedback" oper- 
ation of the network. This is equivalent to computing 
Eq. (60) with nonzero mean, /i. In the specific case of ori- 
entation selectivity, with n = 2, the final distribution of 
tuning strength is then obtained as 

0 2 tecr 2 )=4 e -^ 2 ^ 2 )/^ 2 )/o(^), (63) 
a z a 1 

where 

/ 0 (z) = - f 1 e zcos(0) d9 
x Jo 

is the modified Bessel function of the first kind and order 
zero. 

Endnotes 

a Feature-specific connectivity, however, appears later 
during deveolpment (Ko et al. 2013), consistent with 
previous reports (Jia et al. 2010; Ko et al. 2011). 

b Note that, as a TT-periodic parameter 9 is now mapped 
to a 2jr -periodic parameter 29, an opposite direction 
here implies an orthogonal orientation in the original 
parameter space. 

c One should also consider the possible effect of 
feature-specific connectivity (Ko et al. 2011, 2013) on this 
behavior, which predicts a preponderance of connections 
between neurons with similar selectivity and hence a 
more tuned input to neurons. 

d Tuning width is unchanged if the baseline activity is 
removed, and half-width at half-height (HWHH) is 
computed from the Gaussian fit to the tuned part, as 
done in (Sharon and Grinvald 2002). Tuning 
amplification would effectively decrease the tuning 
width, however, if the baseline is also taken into 
consideratio (as in Wilson et al. (2012)). 

e Note that we do not exclude the presence of nonlinear 
mechanisms in the biological cortex and their 
contribution to orientation selectivity. The nonlinearity 
of the neuronal transfer function, as well as other 
nonlinear mechanisms like nonlinear dendritic 



amplification (Jia et al. 2010; Lavzin et al. 2012; Lee et al. 
2012a) or synchronization of thalamic inputs (Bruno and 
Sakmann 2006; Stanley et al. 2012), may contribute in 
addition to obtain higher selectivity. 

f We verified however, if the cosine fit imposes a general 
constraint on the tuning curves and changes the OSI*. To 
check this, we use an alternative way to obtain /"p re f and 
r ort h from the smoothened tuning curves, namely by 
linear interpolation between the data points. The result 
of this is compared with the result of the cosine fit in 
Figure 3B, and does not change the conclusions 
qualitatively. 

g If there is no minimum boundary condition, as it is the 
case for our neuron model, this hyperpolarization can 
grow beyond all bounds. Note that this would not be the 
case if one uses another neuron model, like 
conductance-base LIF neurons. The treatment of the 
problem in those cases might, therefore, be different. 
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